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ABSTRACT 


A  STUDY  OF  MOVEMENT  DETECTION  IN  FUNCTIONAL  ECHO-PLANAR 
IMAGING  OF  THE  BRAIN.  Ernst  C.  Hansch,  Gregory  McCarthy. 
Section  of  Neurosurgery,  Department  of  Surgery,  Yale 
University,  School  of  Medicine,  New  Haven,  CT. 

Functional  echo-planar  MRI  imaging  of  the  brain  is  limited  by 
patient  movement  between  images.  This  study  examined  the 
ability  to  detect  and  quantify  patient  movement  by  analysis 
of  echo-planar  images  of  the  head  themselves,  without  external 
fiducial  markings  or  outside  references.  A  center  of  mass 
algorithm,  a  comparison  of  the  variance  of  ratios  of  pixel 
intensities,  and  an  area  comparison  algorithm  were  all 
evaluated  initially  with  a  computer  model  of  movement,  then 
with  a  set  of  actual  perturbed  echo-planar  images,  and  finally 
with  echo-planar  images  of  a  phantom  model.  Center  of  mass 
changed  with  movement,  but  was  limited  in  its  application  to 
quantitation  of  linear  translation  in  the  imaging  plane  (R2  of 
0.9972).  The  ratios  of  variance  algorithm  detected  and 
quantitated  linear  translation  and  rotation  in  the  imaging 
plane  (R2  0.9991  and  0.9865),  but  was  slow  to  apply  and  unable 
to  quantitate  movement  through  the  imaging  plane.  The  area 
algorithm,  while  able  to  detect  movement  through  the  imaging 
plane  (average  area  change  25 . 49%) , could  not  quantitate  the 
movement.  Movement  is  a  serious  problem  in  functional  echo- 
planar  imaging  that  is  not  easily  resolved  by  techniques 
involving  post-processing  of  echo-planar  image  data. 
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INTRODUCTION 

By  observing  individuals  with  brain  injuries  or  lesions, 
investigators  over  the  years  came  to  the  conclusion  that 
specific  functions  are  controlled  by  certain  discreet  areas  of 
the  brain.  As  long  ago  as  the  time  of  Hippocrates, 
observations  were  made  that  hemiplegia  often  resulted  from 
brain  injury  to  the  contralateral  hemisphere1.  Case  histories, 
such  as  the  report  of  a  mine  worker  surviving  passage  of  an 
iron  rod  through  his  frontal  lobe,  associated  higher  human 
function  with  certain  lesions  of  the  brain2.  Over  the  years 
these  early  observations  were  expanded-on  with  studies  based 
on  electrical  stimulation  of  the  brain.  Fritsch  and  Hitzig 
were  able  to  accurately  localize  motor  activity  in  the  dog 
through  such  experiments1.  The  work  of  Broca  and  Wernike 
allowed  for  the  localization  of  speech  activity  in  the  brain3. 
Further  mapping  of  the  cortical  surface  was  performed  by 
Brodmann4  and  Penfield5.  Early  surgical  procedures  were  often 
performed  with  the  patient  awake.  By  stimulating  those  areas 
of  the  brain  being  operated  on,  the  surgeon  could  develop  a 
better  understanding  of  the  impact  of  any  planned  surgical 
resections. 

Localization  of  specific  functional  areas  plays  a  major 
role  in  the  clinical  neurosciences.  While  higher  human 
functions  remain  poorly  localized,  a  good  understanding  has 
developed  of  the  distribution  of  motor  and  sensory 
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modalities6.  This  information  is  frequently  applied  in 
neurological  surgery,  particularly  in  cases  requiring 
resection  of  brain  matter.  Tumor  resections  are  carefully 
planned  to  avoid  resection  of  areas  of  the  motor  strip, 
occipital  cortex,  or  speech  centers.  Epilepsy  surgery  often 
involves  the  resection  of  the  amygdala,  hippocampus  and 
temporal  lobe  tip7.  Previous  localization  of  the  dominant 
hemisphere  and  the  speech  center  is  critical  in  these 
instances,  as  the  resection  of  these  sensitive  areas  can  lead 
to  aphasia,  or  short-term  memory  loss.  The  methods  for  such 
localization  have  thus  far  been  complex  and  largely 
interventional . 

Both  functional  and  dysfunctional  activity  must  be  pin¬ 
pointed  in  epilepsy  surgery.  Localization  of  the  epileptogenic 
focus,  as  well  as  the  elucidation  of  the  functional  activity 
of  the  surrounding  brain  is  necessary.  The  site  of 
epileptogenic  activity  is  often  initially  obtained  by 
electroencephalography.  Dependent  on  the  detection  of  electric 
brain  activity,  this  method  is  one  of  the  oldest  and  most 
widely-used  methods  of  functional  localization8.  More  precise 
localization  of  epileptogenicity  in  the  area  of  interest  can 
be  obtained  through  subdural  cortical  mapping9.  In  this 
procedure,  an  electrode  grid  is  surgically  placed  on  top  of 
the  suspected  epileptogenic  cortex.  The  activity  of  the 
cortical  area  is  thus  monitored  over  a  period  of  time  to 
further  localize  the  epileptogenic  focus.  The  determination  of 
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the  dominant  side  of  the  brain  is  determined  by  the  WADA  test, 
in  which  Sodium  Amytal  is  injected  into  the  left  and  right 
brain  circulation  respectively10.  This  transiently  paralyzes 
the  injected  side,  allowing  observation  of  the  patient's 
resultant  behavior.  Memory  and  speech  changes  thus  implicate 
the  paralyzed  side  as  the  dominant  one.  More  detailed 
localization  of  speech  and  memory  may  be  attained  at  the  time 
of  surgery  by  means  of  electrical  stimulation  of  the  planned 
area  of  resection  with  the  awake  patient11  . 

Other  experimental  techniques  of  interventional 
functional  imaging  are  also  being  evaluated.  Optical  imaging 
methods  based  on  the  detection  of  activity  related  changes  in 
the  optical  reflectance  properties  of  the  cortex  are  under 
investigation.  These  methods  has  thus  far  been  successfully 
applied  in  both  anesthetized12,13  and  awake  animals14. 

Given  the  importance  of  functional  localization  of 
cerebral  activity,  it  is  clear  that  non-invasive  techniques 
for  determination  of  activity  would  be  particularly  useful.  A 
number  of  imaging  modalities  have  proven  to  be  of  benefit  in 
this  regard.  These  include  positron  emission  tomography  (PET) , 
single-photon  emission  computed  tomography  (SPECT) ,  and 
nuclear  magnetic  resonance  techniques  (both  magnetic  resonance 
spectroscopy  [MRS]  and  magnetic  resonance  imaging  [MRI ] ) . 


PET  and  SPECT 
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Positron  Emission  Tomography  (PET)  is  based  on  the 
knowledge  of  localized  metabolism  and  blood  flow  to  identify 
regions  of  activity15,16.  A  radioactive  tracer  is  introduced,  and 
its  distribution  within  the  brain  monitored  by  its  decay 
through  emission  of  positrons.  The  tracer  is  produced  by  the 
incorporation  of  an  isotope  (e.g.  150,  nC,  13N,  18F)  into  a  compound. 
Metabolic  activity  dependent  on  the  compound  may  thus  be 
monitored.  For  example,  a  bolus  of  intravenous ,  150-labelled 
water  is  injected  into  the  study  subject,  and  any  activity- 
related  changes  in  cerebral  blood  flow  is  measured.  PET 
imaging  has  seen  relatively  widespread  clinical  use,  being 
used,  amongst  other  fields,  in  the  study  of  psychiatric 
disorders17,  epilepsy18,  Parkinson's  disease19  and  brain  tumors20. 
It  has  also  been  widely  used  in  the  academic  neurosciences. 
Using  this  technique,  Petersen  et.al.  found  activation  in 
certain  areas  of  the  left  medial  extrastriate  visual  cortex 
upon  visual  representation  of  words21,  and  found  evidence  that 
the  various  codes  of  word  processing  are  localized  in  widely 
varying  areas  of  the  brain22.  The  same  group  has  also  done  work 
in  localizing  activity  in  the  visual  cortex.23  More  complex 
frontal  lobe  dependent  processes  have  also  been  studied  by 
these  techniques.  In  a  study  of  word  finding,  Frith  et.  al. 
discovered  an  increase  in  activity  in  left  dorsolateral  pre¬ 
frontal  activity  with  tasks  stimulating  intrinsic  word 
generation  24 . 

Despite  these  success,  PET  studies  suffer  from  a  number 
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of  limitations.  PET  scanners  are  large  and  expensive,  limiting 
their  availability.  Because  of  their  dependence  on  radioactive 
tracer  materials,  they  must  generally  be  located  close  to 
cyclotrons  where  such  tracers  are  produced.  Furthermore,  the 
exposure  to  radioactive  tracers  limits  the  use  of  PET  studies 
in  any  one  individual.  Many  repeat  studies  cannot  be  performed 
due  to  the  resultant  radiation  exposure.  PET  images  are  of 
relatively  low  resolution  (5.0mm  X  5.0mm  X  5.0mm  voxels), 
limiting  the  accuracy  of  functional  localization.  Also, 
because  of  the  low  signal  to  noise  ratio  of  PET  images,  the 
results  of  several  subjects  must  often  be  pooled  to  obtain 
satisfactory  results. 

Like  PET  imaging,  Single-Photon  Emission  Computed 
Tomography  (SPECT)  is  based  on  the  estimation  of  activity  by 
the  analysis  of  radioactive  tracers25,26.  The  isotopes  used  in 
SPECT  imaging,  however,  emit  photons  rather  than  positrons  and 
are  therefore  easier  to  detect.  Because  of  the  lipophilic 
qualities  of  most  SPECT  tracers,  and  their  resultant  ability 
to  cross  the  blood-brain  barrier,  SPECT  is  frequently  used  for 
perfusion  imaging27.  Clinical  application  have  included  studies 
of  strokes28,  blood-brain  permeability29,  and  epilepsy30.  While 
SPECT  imaging  monitors  perfusion  directly,  it  can  also  be  used 
to  provide  information  regarding  metabolic  activity  of  the 
brain31 . 

Because  of  the  relative  ease  of  detection  of  SPECT  tracer 
decay,  SPECT  imaging  is  cheaper  and  simpler  to  use  than  PET 
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imaging.  Furthermore,  due  to  the  longer  half-life  of  the 
isotopes  used,  it  is  also  not  as  strictly  dependent  on 
proximity  to  a  cyclotron.  However,  this  imaging  modality 
nevertheless  suffers  from  a  number  of  the  shortcomings  of  PET 
imaging.  The  resolution  obtained  is  actually  lower  than  that 
of  PET.  SPECT  also  exposes  the  subject  to  radiation,  thereby 
limiting  its  application  in  a  single  subject. 


Nuclear  Magnetic  Resonance  Technique 

Recently,  magnetic  resonance  imaging  techniques  have  been 
developed  which  hold  out  the  promise  of  providing  simple,  non- 
invasive  methods  of  functional  imaging.  Developed  initially  in 
197332,  the  first  clinically  useful  MRI  machines  began 
appearing  in  1983.  While  allowing  for  excellent  anatomical 
soft-tissue  differentiation,  MRI  techniques  have  long  suffered 
from  the  drawback  of  relatively  long  image  acquisition  times, 
limiting  their  use  to  those  areas  of  the  human  body  where 
little  or  no  motion  occurs.  This  is  the  case  with  the  spin- 
echo  (SE)  sequence,  the  most  widely  used,  but  slow,  imaging 
sequence.  Improvements  in  speed  from  many  minutes  to  a  few 
seconds  came  in  the  form  of  the  gradient-echo  techniques, 
first  used  in  the  mid  80's33,34.  Echo-planar  imaging  (EPI)  ,  the 
technique  that  has  turned  out  to  have  revolutionary  effects  on 
functional  magnetic  resonance  imaging,  was  initially  described 
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in  1977 ,  by  Peter  Mansfield35.  Though  only  recently  applied, 
this  technique  underwent  further  laboratory  development  in  the 
80's  36’37,  and  is  now  beginning  to  undergo  clinical  use. 

Magnetic  resonance  imaging  takes  advantage  of  the 
property  of  angular  momentum  possessed  by  elements  with  an  odd 
atomic  number38.  Due  to  its  abundance  in  the  human  body,  the 
hydrogen  nucleus  ^H)  is  generally  used  in  medical  MR  I 
imaging.  The  magnetic  field  associated  with  the  spinning 
hydrogen  nucleus  is  referred  to  as  the  magnetic  dipole  moment. 
By  application  of  a  large  external  magnetic  field  to  the 
subject  being  imaged,  an  alignment  of  the  magnetic  dipole 
moments  of  the  subject's  hydrogen  nuclei  occurs.  The  resultant 
magnetic  dipole  moment  vector  may  be  thought  of  as  the 
subject's  magnetization,  and  is  referred  to  as  the  B0  vector. 
By  application  of  the  external  magnetic  field,  the  magnetic 
dipole  moments  begin  to  precess  at  a  frequency  given  by  the 
Larmor  Equation: 


where  y  is  the  constant  gyromagnetic  ratio  for  hydrogen,  and 
G0  the  externally  applied  magnetic  field.  The  frequency  of  the 
procession  is  directly  proportional  to  G0.  By  further 
application  of  a  series  of  radiofrequency  pulses,  with  their 
magnetic  field  component  perpendicular  to  the  direction  of  the 
external  magnetic  field,  a  spiraling  of  the  B0  vector, from  the 
y-axis,  towards  the  x-z  plane  is  attained  (Figure  1) . 
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Figure  1 
Motion  of  the 
B0  field 
towards  the 
x-z  plane, 
while 

precessing 
around  the  y 
axis. 


To  achieve  this,  the  RF  pulses  must  be  applied  at  the 
resonant  Larmor  frequency  of  the  hydrogen  nucleus.  As  the  B0 
field  spirals  towards  the  x-z  plane,  it  may  be  detected  as  a 
magnetic  resonance  signal  in  the  form  its  component  Bj  in  the 
x-z  plane.  Finally,  upon  the  termination  of  the  RF  pulse,  a 
relaxation  occurs  characterized  by  the  gradual  return  of  the 
B0  field  along  the  y-axis.  The  rates  of  relaxation  are  given 
by  the  time  constants  Tlf  and  T2,  which  correspond  respectively 
to  the  rate  of  growth  of  the  B0,  vector  and  the  rate  of 
shrinkage  of  the  Bt  component.  The  magnetic  resonance  signal 
decays  with  the  decay  of  B{.  Through  the  use  of  gradient 
coils,  a  linear  relationship  may  be  established  between  the 
applied  external  magnetic  field  and  position.  By  pulsing  these 
coils  on  the  x,y  and  z  axes,  the  magnetic  resonance  signal  may 
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be  spacially  encoded,  and  an  image  generated  by  the 


application  of  a  2  dimensional  Fourier  transform  : 

Sin,  t)  =  J [m{x,  y)  *V'dxdy 

Mix.y)  =J Js[n,  t)  e^^e'^^dndt 


Where  M(x,y)  represents  the  magnetization  of  the  pixel  at 
point  x,y,  and  S(n,t)  represents  the  MR I  signal. 

Magnetic  resonance  spectroscopy  (MRS)  uses  the  principle 
of  magnetic  resonance  to  measure  the  quantity  and  location  of 
specific  metabolites39.  The  resonance  frequency  of  a  certain 
molecule  (e.g.  1H,  3lP,  13C)  will  vary  with  different  molecular 
locations,  even  in  the  same  G0  field.  The  location  of 
metabolites  such  as  lactate,  glycogen  and  ATP  may  therefore  be 
mapped  by  the  identification  of  sources  of  1H,  3lP  and  13C 
respectively.  Areas  of  ischemia  have  been  localized  and 
measured  with  MRS  by  identification  of  lactate  distribution40,41. 
Studies  of  glioma  metabolism  have  been  performed  by 
characterization  of  glioma  glucose  and  lactose  processing42,43  . 
Furthermore,  preliminary  studies  have  held  out  the  promise  of 
identification  of  areas  of  activation  by  monitoring  of  glucose 
metabolic  rate  by  MRS44. 

The  most  dramatic  breakthrough  in  the  field  of  functional 
MRI  imaging,  however,  has  come  in  the  form  of  echo-planar 
imaging.  Echo  planar  imaging  techniques  are  based  on  the  same 
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principles  as  regular  MRI,  with  a  few  minor  variations.  Unlike 
conventional  MRI  imaging  sequences,  echo  planar  techniques 
require  only  a  single  RF  pulse,  after  which  the  image  data  is 
acquired45.  This  decreases  acquisition  time  to  the  range  of  4  0- 
128  ms.  The  equipment  required  for  EPI  is  similar  to  that 
needed  for  conventional  MRI,  with  a  few  exceptions.  Special  RF 
coils  are  needed  for  higher  magnetic  field  strengths,  shielded 
gradient  coils  are  required  for  rapid  switching  of  the 
external  magnetic  gradient,  and  shim  coils  may  be  needed  for 
minimizing  magnetic  field  variations. 

The  short  imaging  time  possible  with  EPI  has  allowed  for 
the  possibility  of  studying  physiological  phenomena  on  the 
same  short  time  scale.  As  in  PET  imaging,  physiological 
changes  in  the  area  of  interest  can  be  used  as  indicators  of 
activity.  Variations  in  neuronal  activity  are  reflected  in 
changes  in  metabolism  and  perfusion  of  the  surrounding  brain 
tissue46.  In  particular,  changes  in  cerebral  blood  flow(rCBF) , 
cerebral  oxygen  metabolic  rate  (rCMR02)  ,  and  cerebral  blood 
volume  (rCBV)  were  found  to  relate  to  the  relaxation  times  T2* 
and  Tt  47 ,  T2*  being  a  decay  constant  related  to  T2,  but 
shortened  by  inhomogeneities  in  the  externally  applied 
magnetic  field.  It  is  given  by: 
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here  1/T2m  and  l/t2  are  signal  decay  rates  due  to  macroscopic 
and  microscopic  external  field  inhomogeneities. 

Functional  MRI  Studies 

A  landmark  study  by  Belliveau  et.  al.  was  amongst  the 
first  to  detect  functional  activity  by  MRI  techniques. 
Belliveau  and  colleages  were  able  to  demonstrate  activity  in 
the  occipital  cortex  of  humans  by  detecting  an  increase  in 
cerebral  blood  volume  during  visual  stimulation48.  Earlier 
animal  studies  had  determined  that  the  administration  of  0.5M 
gadolinium  diethylenetriaminepentaacetic  acid  [Gd (DTPA) 2] ,  a 
paramagnetic  contrast  agent,  as  a  blood  tracer,  could  lead  to 
creation  of  maps  of  cerebral  blood  volume49,50,51  .  The  presence  of 
paramagnetic  contrast  agent  in  the  capillaries  of  the  brain 
induces  local  microscopic  magnetic  field  gradients  which  in 
turn  reduce  T/,  T2*,  and  T2.  These  changes  could  then  in  turn 
be  detected  by  an  EPI  imaging  sequence  sensitive  to  such 
variations.  A  bolus  of  Gd(DTPA)  injected  intravenously  could 
mix  with  all  the  circulating  blood  in  the  heart,  and  arrive  at 
the  brain  approximately  15  seconds  after  injection.  An 
appropriately  sensitive  EPI  imaging  sequence  could  then  detect 
the  presence  of  the  tracer  in  the  brain.  By  digital 
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subtraction  of  CBV  maps  of  the  resting  state  from  CBV  maps 
during  times  of  task  activation,  areas  of  enhancement  could  be 
obtained.  By  subjecting  subjects  to  photic  stimulation, 
Belliveau  found  an  increase  in  cerebral  blood  flow  of  32 
percent  in  areas  of  the  occipital  cortex  during  times  of 
stimulation48. 

This  technique,  while  allowing  higher  resolution  imaging 
than  available  with  PET  scanners  without  exposure  to 
radiation,  still  suffered  from  dependence  on  Gd(DTPA)  contrast 
administration.  Due  to  the  toxicity  of  the  contrast,  and  the 
delay  in  delivery  to  the  brain,  the  number  of  studies  possible 
on  any  one  individual  were  limited.  Later  studies  were  able  to 
successfully  use  the  properties  of  blood  as  a  physiological 
contrast  material.  It  was  noted  that  hemoglobin  contained 
materials  with  magnetic  properties  -  deoxyhemoglobin  with 
paramagnetic  iron,  and  oxyhemoglobin  with  diamagnetic  oxygen- 
bound  iron52.  The  difference  in  magnetic  properties  between 
deoxyhemoglobin  and  oxyhemoglobin  can  lead  to  microscopic 
variations  in  the  magnetic  field  gradient,  which  may  be 
detected  as  small  changes  in  T2*.  By  using  an  EPI  sequence  with 
a  7T  magnet  on  a  rat,  Ogawa  et.al.  found  that  variations  in 
blood  oxygenation  could  therefore  be  detected  by  changes  in 
magnetic  properties  of  hemoglobin  53,54  .  Further  animal  studies 
by  Turner  et.al.  using  a  long-TE  EPI  sequence  sensitive  to  T2* 
found  variations  in  image  intensity  with  changes  in  hemoglobin 
concentration  in  the  blood55.  This  demonstrated  the  ability  of 
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EPI  to  follow  transient  oxygen  changes  in  the  blood. 

Further  human  studies  bore  out  the  potential  of  the  new 
hemoglobin-based  technique.  Using  this  technique,  Bandettini 
et.  al.  found  a  signal  increase  of  4.3%  in  the  primary  motor 
and  sensory  cortices  during  a  finger  moving  task56.  Kwong  et. 
al.  were  able  to  use  both  gradient  echo  (GE)  and  spin-echo 
inversion  recovery  (IR)  imaging  sequences  to  observe  changes 
in  the  primary  visual  cortex  of  subjects  during  photic 
stimulation57.  By  using  a  4T  magnet  with  a  gradient-echo 
imaging  sequence,  Ogawa  et.  al.,  were  able  to  more  easily 
detect  signal  changes,  and  thereby  decrease  their  image  voxel 
size.  They  found  easily-detectable  transient  signal  changes 
confined  to  the  gray  matter  the  primary  visual  cortex58. 
Blamire  et.al.  used  a  2.1  T  magnet  to  study  activation  in  the 
visual  cortex,  finding  an  increase  in  signal  of  9.7%,  but 
delayed  by  3.5  -  5.0  seconds  after  the  onset  of  photic 
stimulation59.  Similar  delays  were  noted  in  studies  by  Kwong57. 
A  study  by  McCarthy  and  colleagues  investigated  the 
distribution  of  language-based  tasks.  By  having  subjects 
generate  and  read  words  in  response  to  words  presented  to  them 
by  the  experimenter,  this  group  was  able  to  note  increased 
signal  in  a  sulcus  anterior  to  the  lateral  sulcus,  Brodman's 
area  47  and  area  1060. 

It  is  clear  that  while  the  technique  of  hemoglobin-based 
functional  MRI  imaging  is  still  in  its  infancy,  it  has  obvious 
advantages  over  other  methods  of  functional  imaging.  It  is 
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capable  of  producing  higher  resolution  images  than  PET 
studies.  The  absence  of  contrast  allows  the  experiments  to  be 
carried  out  on  the  same  subject  numerous  times,  with  ease  and 
safety.  Due  to  the  improved  S/N  signal  ratio  than  available 
with  PET  imaging,  only  one  subject  is  required  for  a  complete 
study61.  The  results  may  therefore  be  more  easily  used  in  the 
context  of  surgical  planning. 

Movement  Artifacts 


The  technique  is  not  without  its  problems,  however.  The 
processing  of  functional  MR  I  data  is  based  on  a  digital,  voxel 
by  voxel  analysis  of  MRI  images  to  determine  the  area  of 
activation.  In  the  course  of  a  typical  experiment,  a  series  of 
images  is  made  under  different  conditions  of  subject 
stimulation.  For  example,  in  a  typical  experiment  as  performed 
by  the  group  at  the  West  Haven  VA,  a  series  of  60  images  is 
taken.  Twenty  are  taken  before  the  initiation  of  the  task, 
twenty  with  the  subject  performing  a  certain  task,  and  twenty 
after  the  completion  of  the  task.  The  resultant  images  are 
then  compared,  and  those  taken  with  the  patient  in  the  resting 
state  digitally  subtracted  from  those  taken  at  the  time  of 
task  activation.  In  particular,  a  t-test  is  applied  to  the 
images  on  a  pixel  by  pixel  basis,  searching  for  significant 
changes  in  pixel  intensity.  A  t-test  map  of  significantly 
different  pixels  for  a  chosen  value  of  T  is  thereby  produced, 


. 
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with  the  areas  of  significant  change  thought  to  correlate  with 
areas  of  activity.  These  maps  may  be  superimposed  on  an 
anatomic  MR I  image  of  the  subject  made  before  the  begin  of  the 
functional  study,  thereby  producing  an  anatomic  map  of 
activity  (see  figure  2). 


Figure  2 

Activation  of  the 
motor  cortex  in  a 
finger-moving  task. 
Activation  results 
are  superimposed  on 
an  anatomical 
image . 


For  the  resultant  image  to  be  of  value  in  assessing  areas 
of  activation,  the  supposition  is  made  that  no  motion  of  the 
imaging  subject  occurs  between  the  images  taken  during  task 
activation  and  of  the  images  at  baseline.  If  motion  does 
occur,  the  resultant  subtracted  image  could  contain 
significant  artifactual  errors.  Restraint  of  the  imaging 
subject,  along  with  instructions  to  remain  still  during  the 
imaging  sequence  may  be  helpful  in  minimizing  movement.  In  a 
typical  experiment  by  the  West  Haven  group,  the  subject's  head 
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is  placed  in  a  vac-pac  surgical  positioning  mattress  restraint 
system,  cushioned  with  styrofoam  to  reduce  motion.  The  air  is 
removed  from  the  vac-pac  to  fix  it  in  position,  and  a  band  is 
placed  across  the  subject's  forehead  to  maintain  the  head  as 
still  as  possible.  Despite  these  precautions,  significant 
motion  may  still  occur.  The  series  of  images  may  take  as  long 
as  an  hour,  during  which  even  the  best-intentioned  subjects 
often  move  gradually.  The  task  itself,  if  it  requires  the 
pressing  of  a  button,  the  opening  and  closing  of  eyes,  or  the 
moving  of  the  jaw  can  induce  the  subject  to  move  against  the 
restrains.  Even  in  the  most  optimistic  scenario,  it  must  be 
remembered  that  a  volunteer  study  subject  is  more  likely  to 
remain  still,  and  be  willing  to  repeat  the  study  if  needed, 
than  is  an  actual  patient  having  the  study  performed  for 
diagnostic  reasons.  As  the  technique  becomes  more  widely 
available,  and  becomes  more  frequently  used  in  a  clinical 
setting,  the  limitations  of  movement  artifacts  will  become 
more  acutely  evident. 

The  types  of  movement  possible  are  easily  described  in 
geometrical  form.  Given  the  familiar  3-dimensional  axis  system 
(see  figure  3),  the  patient's  motion  may  be  seen  as  a 
combination  of  motions  along  and  around  the  3  orthogonal  axes. 
Defining  the  imaging  plane  as  the  plane  of  reference,  one 
could  describe  patient  motion  as  being  either  in  the  plane  of 
imaging,  through  the  plane  of  imaging,  or  a  combination  of  the 


two. 
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Figure  3 
Orthogonal 


axis  system.  The 
imaging  plane  is 
defined  to  lie  in 
the  x-y  plane. 


For  example,  in  a  patient  lying  in  an  MR I  magnet  having 
axial  images  taken,  rotation  of  his  head  side-to-side  or 
shifting  of  his  head  side-to-side  would  produce  a  movement  in 
the  plane  of  imaging.  Similarly,  a  nodding  motion  would 
produce  motion  through  the  plane  -  a  combination  of  rotation 
and  translation  through  plane. 

The  effect  of  patient  movement  may  be  seen  in  the 
analysis  of  t-maps,  as  described  above.  In  moving,  areas  of 
activation  are  shifted  in  the  direction  of  motion.  Pixels  at 
the  borders  of  contours  or  on  the  edge  of  the  activated  area 
therefore  begin  to  change  their  intensity  to  correspond  with 
the  intensity  changes  at  the  contour  edge.  A  rim  of  intensity 
changes  of  significantly  different  t-values  therefore  appears 
around  the  image.  The  resultant  "edge  effect"  as  seen 
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illustrated  in  Figure  4,  is  often  a  sign  of  movement  artifacts 
in  the  images.  Such  artifact  makes  identification  of  actual 
activation  difficult. 


Figure  4 

Edge  effect  seen 
superimposed  on  an 
anatomical  image. 
Intensity  changes 
at  contour  borders 
are  enhanced  by 
subject  movement, 
causing 
significant 
artifact . 


A  quantification  of  the  amount  of  error  as  a  function  of 
movement  in  the  plane  may  be  seen  in  Figure  5,  where  error  is 
given  as  a  percentage  of  the  area  of  activation.  By  studying 
these  graphs,  it  becomes  evident  that  the  amount  of  error 
depends  on  a  number  of  factors,  including  the  size  of  the 
area  of  interest  and,  in  the  case  of  rotation,  the  distance  of 
the  area  of  interest  from  the  point  of  rotation.  In  general, 
movement  results  in  a  linearly  increasing  error.  As  would  be 
expected,  because  of  their  small  size,  small  areas  of 
activation  are  subject  to  more  rapid  increases  in  error  with 
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movement  (Fig  5.  a).  In  the  case  of  rotation,  the  rate  of 
increase  depends  on  the  distance  of  the  area  of  interest  from 
the  point  of  rotation,  with  faster  increase  in  error  at 
greater  distances  (Fig  5.b)  .  If  we  view  the  error  due  to 
motion  in  the  plane  as  a  vector,  with  components  Ex,  Ey,  Er 
corresponding  to  translation  along  x,  translation  along  y  , 
and  rotation  in  plane  respectively,  we  can  calculate  the  total 
error  by: 

E=sjEl+E*+E2r 

In  figure  5.c  we  may  see  the  total  amount  of  movement  for  a  20 
X  15  pixel  elliptically  shaped  area  of  activation  at  50  pixels 
from  the  point  of  rotation.  As  may  be  seen,  the  total  error 
increases  linearly  with  the  amount  of  movement.  Movement  of 
the  patient  by  7  pixels  along  both  the  x  and  y  axes,  combined 
with  rotation  of  7  degrees  around  the  z  axis  results  in  an 
error  as  large  as  the  size  of  the  area  of  activation.  Even 
with  motion  of  1  pixel  along  each  axis  and  1  degree  of 
rotation,  the  total  error  is  approximately  15%  of  the  total 
activated  area. 

Clearly,  therefore,  a  method  of  dealing  with  movement  is 
necessary  for  the  proper  analysis  of  functional  MRI  data.  The 
physical  restraint  of  the  patient  has  already  been  mentioned 
as  a  possible  solution.  While  this  approach  may  be  effective 
in  most  situations,  it  does  not  eliminate  all  movement. 
Analysis  of  the  image  data  for  movement,  with  the  intent  to 
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detect  and  possibly  correct  for  the  movement  would  certainly 
be  extremely  helpful.  The  aim  of  such  analysis  would  be  to 
detect  motion  on  a  scale  of  a  few  millimeters.  For  an  image 
displayed  on  a  computer  screen,  this  correlates  to  detection 
on  the  order  of  approximately  2  pixels.  Because  the  patient 
is  physically  restrained,  he  or  she  is  unable  to  move  more 
than  a  few  millimeters  or  degrees.  Furthermore,  motions  of 
more  than  a  few  millimeters  may  be  detected  by  visual 
inspection  of  the  resultant  images,  and  the  data  set 
appropriately  dealt  with.  While  this  amount  of  motion  may 
appear  small,  as  noted  above  it  may  result  in  significant 
error. 

Little  work  has  been  done  in  this  regard.  The  work  on 
magnetic  resonance  image  motion  artifact  correction  has  thus 
far  been  aimed  largely  at  correcting  for  patient  motion  during 
the  scanning  acquisition  period.  As  conventional  spin-echo  MR  I 
acquisition  times  are  on  the  order  of  minutes,  this  problem 
has  been  a  source  of  image  degradation  in  the  past,  and  work 
has  been  done  in  reducing  the  amount  of  ghosting  and  blurring 
of  such  images62.  Due  to  the  rapid  acquisition  time  of  EPI, 
however,  this  work  does  not  apply  to  the  problem  of  patient 
motion  between  images.  Other  than  restraint  of  the  study 
subject,  no  other  systematic  approach  to  functional  MRI 
imaging  movement  control  or  detection  has  been  reported.  It  is 
the  purpose  of  this  work  to  study  this  problem.  In  particular, 
an  attempt  will  be  made  to  detect  movement  by  application  of 
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post-processing  techniques  on  the  echo-planar  images 
themselves,  independent  of  any  external  fiducials  or  external 
references . 
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Figure  5 .a 

Error  as  a 
function  of 
activation 
area  size 
for  constant 
size  shift. 


Figure  5 . b 

Error  as  a 
function  of 
distance  of 
area  of 
activation 
from  point 
of  rotation 


Figure  5.c 

Total  error 
of  a  20  X  15 
pixel  area  of 
activation, 
and 

components . 
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STATEMENT  OF  PURPOSE: 


The  purpose  of  this  study  was  to  evaluate  the  feasibility 
of  movement  detection  in  functional  echo-planar  imaging  by  the 
development  and  testing  of  computer  algorithms  for  that 
purpose.  In  particular,  an  attempt  was  made  to  limit  the 
movement  detection  techniques  to  those  involving  processing  of 
the  echo-planar  images  themselves,  without  reliance  on  any 
further  reference  frames  or  fiduciary  markings.  Non-reliance 
on  fiduciary  markings  would  allow  any  developed  techniques  to 
be  applicable  retrospectively. 


Page  24 


METHODS  AND  THEORY 


The  problem  of  movement  was  studied  by  development  of 
post-processing  computer  algorithms.  These  were  applied  and 
tested  in  three  levels  of  increasing  realism  -  a  computer 
model  of  movement,  an  analysis  of  perturbed  actual  echo-planar 
images,  and  analysis  of  echo-planar  images  of  a  phantom  model. 
All  programming  code  was  written  in  C++  using  the  Zortech© 
C++  compiler  with  Flash  Graphics©  on  an  IBM  PC.  All 
programming  was  performed  by  the  author,  with  discussion  and 
evaluation  with  Dr.  Gregory  McCarthy,  Anthony  Adrignollo,  and 
Frances  Favorini.  Conversion  of  the  EPI  data  from  its  initial 
format  to  a  256  X  256  floating  point  array  format  as  read  by 
the  author's  programs  was  performed  by  the  MR  software  package 
developed  for  the  neuropsychology  laboratory  at  the  West  Haven 
VA  hospital  by  Anthony  Adrignollo. 

I.  Image  Processing  Techniques 

Upon  generation  and  display  of  the  images,  a  number  of 
image  processing  techniques  were  programmed  to  attempt  to 
detect  image  motion.  The  motion  itself  was  separated  into  two 
components.  Any  changes  in  the  image  were  analyzed  from  the 
perspective  of  motion  in  the  plane  of  imaging,  and  motion 
outside  of  the  plane  of  imaging.  Defining  the  reference 
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system  around  the  MR  I  image,  the  imaging  plane  was  designated 
as  the  z=0  (x-y)  plane,  with  its  origin  at  the  center  of  the 
image  (128,128),  and  the  z-axis  coming  out  of  the  screen 
towards  the  viewer  (Figure  3) .  Motion  in  the  plane  of  imaging 
was  thus  analyzed  as  translation  along  the  X  and  Y  axes,  and 
rotation  around  the  Z  axis.  Similarly,  motion  through  the 
imaging  plane  was  seen  as  translation  along  the  Z  axis,  and 
rotation  around  the  X  and  Y  axes.  This  convention  was  followed 
throughout . 

Image  processing  techniques  were  implemented  in  analysis 
of  the  data.  These  were  performed  in  the  form  of  convolutions 
by  template  operations  for  the  purpose  of  noise  removal,  or 
edge  detection  (low  pass  filters  and  Laplacian  filters) .  Given 
an  image  I(X,Y)  of  size  N,M,  the  convolution  of  the  image  by 
a  template  T(x,y)  of  size  (n  X  m)  was  calculated  by63. 


The  standard  3X3  low  pass  filter: 


111 

111 

111 


was  used  for  noise  removal,  while  the  Laplacian  templates 
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were  used  for  the  edge  detection  techniques. 

The  algorithms  of  motion  detection  were  developed  and 
evaluated  from  the  standpoint  of  motion  detection  in,  and  out 
of  the  plane  of  imaging.  Three  general  methods  were 
implemented:  a)  comparison  of  image  areas;  b) comparison  of 
image  center  of  mass;  c)  comparison  of  the  ratio  of  pixel 
intensity.  The  theory  behind  these  will  be  discussed  in  the 
following  sections. 

I.  a.  Area  Comparison 

The  tested  object  is  imaged  through  an  imaging  plane, 
producing  a  cross-sectional  view.  Through  movement,  the  plane 
of  intersection  is  changed,  producing  different  cross- 
sectional  views  of  the  object.  By  monitoring  the  area  of  the 
resultant  cross-sectional  image,  an  indication  of  motion  could 
be  obtained.  In  particular,  in  the  case  of  motion  in  the 
imaging  plane,  no  change  in  area  is  expected  to  occur.  In 
moving  through  the  imaging  plane,  however,  a  change  in  area 
would  be  expected,  commensurate  with  the  degree  of  motion. 

The  procedure  for  calculation  of  area  was  designed  to  be 
independent  of  the  shape  of  the  object  being  processed. 
Irregularly  shaped  images  could  have  areas  calculated  as  could 
mathematically  well  defined  outlines.  The  procedure  was 
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implemented  by  the  application  of  edge  detection  and  filling 
algorithms.  The  initial  image  would  be  processed  by  a  low  pass 
filter  to  remove  noise  and  thresholded  to  accentuate  the 
edges.  A  Laplacian  operator  as  described  above  would  then  be 
applied  to  delineate  the  edges  if  necessary.  With  the  edges 
of  the  image  thus  outlined,  the  area  could  be  calculated  by 
the  application  of  a  flood  filling  algorithm63. 


I.b.  Center  of  Mass 

Given  a  group  of  particles,  each  with  a  mass  irij,  and 
vector  position  ru  the  total  mass  of  the  system  is  given  by 


in  ■ 


and  the  center  of  mass  by 


m ir 


i 


The  center  of  mass  is  a  particularly  valuable  concept  in 
physics,  for  the  whole  system  moves  like  one  particle  located 
at  the  center  of  mass64.  Motion  of  the  whole  system  may  be 
monitored  by  following  the  motion  of  the  center  of  mass  - 
independent  of  the  motion  of  any  individual  particle  in  the 
system.  Any  external  forces  acting  on  the  system  act  as  one 
force  acting  on  the  center  of  mass. 

An  attempt  was  made  to  extend  this  concept  to  the 
problem  of  motion  detection  in  images.  The  image  being  studied 
is  stored  in  the  form  of  an  256  X  256  array  of  pixels,  which 
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may  be  thought  of  as  a  system  of  particles.  The  intensity  of 
the  individual  particles  may  be  thought  of  as  their  masses. 
Similarly,  the  position  of  the  pixels  may  be  calculated  as 
their  location  relative  to  the  center  of  the  array  (128,128). 
The  resultant  position  would  be  given  by 

X'=  ( Xarray~ 1 2  ®  ) 

y/=(yarray-128) 

The  center  of  mass  of  the  image  was  thus  calculated.  Changes 
in  the  center  of  mass  were  studied  as  possible  indications  of 
motion. 

Strictly  speaking,  the  sum  of  the  intensities  of  all  the 
pixels  in  the  two  images  would  have  to  be  the  same  for  the 
center  of  mass  approach  to  have  absolute  predictive  value. 
While  this  is  possible  in  an  ideal  scenario  of  a  model,  it  is 
very  unlikely  in  actual  MR  I  images.  Small  variations  in 
transmission  and  processing  of  images  would  result  in  small 
intensity  variations  in  two  images  taken  of  identical  areas. 
However,  the  assumption  may  be  made  that  the  variations  will 
be  small,  thereby  allowing  a  good  estimate  of  movement  with 
the  center  of  mass. 


J.c.  Intensity  Ratio 
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This  technique  is  based  on  a  method  developed  by  Woods 
et.  al.  for  alignment  of  3-D  PET  and  MR  I  data  sets65.  In  their 
study,  Woods  and  colleagues  attempted  to  align  two  3-D  image 
sets  by  comparison  on  a  pixel-by-pixel  basis.  The  algorithm 
calculates  the  ratio  of  intensity  of  corresponding  pixels  in 
the  data  sets.  If  the  intensity  of  pixel  i  in  image  A  is 
designated  by  A;,  and  by  in  image  B,  then  the  ratio  r  of  the 
intensities  is  given  by 


The  mean  of  the  ratios  of  N  pixels  may  be  calculated  by  the 
familiar  equation: 


Similarly,  the  standard  deviation  of  the  pixels  is  given  by63 


0  = 


N  N- 


By  dividing  the  standard  deviation  of  the  ratios  by  the  mean 
of  the  ratios,  one  obtains  a  measure  of  the  variance  v: 


In  a  perfectly  aligned  image,  as  a  approaches  0,  so  does  the 
value  of  the  variance  v.  Larger  values  of  v  imply  imperfect 
alignment. 
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The  alignment  of  two  3-D  data  sets  may  be  viewed  as  a 
variant  of  the  problem  of  motion  detection  in  MR I  imaging.  An 
attempt  could  be  made  to  align  the  two  images  in  plane.  As 
alignment  is  attempted,  the  amount  of  shift  needed  to  align 
the  images  is  monitored.  If  the  images  are  found  to  be  best 
aligned  without  any  shift,  then  it  may  be  assumed  that  no 
motion  has  occurred  between  the  two  images.  Conversely,  if 
shifting  of  the  images  is  necessary  for  alignment,  then  it  may 
be  assumed  that  the  images  are  separated  by  the  amount  of 
shift  needed  for  alignment. 

With  this  in  mind,  the  algorithm  was  implemented  as 
follows : 

l.One  image  is  designated  as  the  original  image,  while 
the  other  is  designated  as  the  shifted  image. 

2.  The  variance  v  is  calculated  for  the  two  image  sets. 

3.  The  image  set  designated  as  the  original  is 
respectively  translated  one  unit  forward  and  backward 
along  the  x  and  y  axes,  as  well  as  rotated  one  angular 
unit  clockwise  and  counterclockwise  about  the  z  axis.  For 
each  such  transformation  (  a  total  of  6  -  4  translations 
and  2  rotations)  ,  a  new  value  of  v  is  calculated  and 
stored. 

4.  The  process  is  repeated  until  the  smallest  v  value  is 
attained  for  the  desired  area  of  investigation.  Any 
transformation  of  the  images  performed  are  stored  in 


memory . 
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5.  The  transformations  needed  to  achieve  alignment  are 
displayed  on  the  screen.  This  is  the  detected  motion. 

Each  of  the  techniques  described  above  was  evaluated  as 
a  method  for  detection  of  movement.  In  the  case  of  the  model 
ellipsoid  image,  noise  was  added  to  the  image  and  then 
filtered  using  a  low  pass  filter.  The  detection  algorithms 
were  then  applied.  In  each  of  these  cases,  movement  in  the 
plane  of  imaging  as  well  as  movement  through  the  plane  of 
imaging  was  modelled.  The  following  scenarios  were  therefore 
tested. 

Center  of  Mass  in  plane  out  of  plane 
Variance  of  Ratios  in  plane  out  of  plane 
Area  in  plane  out  of  plane 


II.  Ellipsoid  Model 

A  computer  model  was  developed  to  study  the  problem  of 
movement.  Routines  were  written  allowing  for  the  creation  and 
manipulation  of  an  ellipsoid  in  3  dimensions.  The  above 
described  image  processing  techniques  and  algorithms  were 
then  applied  for  the  detection  of  movement  in  the  model.  The 
following  characteristics  were  sought  in  the  model: 


1.  The  modelled  object  was  to  have  an  easily  producible 


’ 


Page  32 


basic  shape,  similar  in  outline  to  the  images  of  the 
head  created  by  MR I . 

2.  The  model  was  to  allow  for  easy  manipulation  of  the 
modeled  object.  Translation  and  rotation  in  and  out  of 
plane  were  to  be  performed. 

3.  The  model  was  to  simulate  the  image  characteristics  of 
an  MRI  image,  with  comparable  image  density,  and  noise 
level . 

4.  The  modeled  object  was  be  stored  and  manipulated  in 
the  same  computer  data  format  as  MRI  images. 

A  3-dimensional  ellipsoid  was  chosen  as  the  object  to  model. 
It  is  well  described  mathematically,  allowing  for  accurate 
reconstruction.  Furthermore,  it's  similar  appearance  to  the 
outline  of  the  human  head  allowed  for  valid  comparison  with 
MRI  images  of  a  human  head  or  phantom. 

The  general  equation  for  an  ellipsoid  with  its  center  at 
xo/Yo/zo/  and  semi  axes  a,b,c  is  described  by: 

(x-x0)2x  (y-y0) 2  x  (z-z0)2  _i 
a2  b2  c 2 

A  cross-sectional  view,  as  would  be  obtained  in  an  MRI,  was 
calculated  by  setting  the  above  equation  equal  to  that  of  a 
plane  of  interest,  given  by: 


Ax+By+Cz+D~Q 


Page  33 

where  A,B,C,D  are  constants.  In  particular,  for  the  study  of 
motion  in  the  plane  of  imaging,  defined  as  z-0 ,  the  above 


equations  give: 


(y-y0)2_1 

b2 


The  ellipse  model  was  designed  to  be  stored  in  the  same  format 
as  an  MRI  image.  The  above  equations  were  therefore  applied  on 
a  point-by-point  basis  to  form  a  256  X  256  array  of  integer 
values,  with  each  point  corresponding  to  a  pixel  intensity 
value  with  the  same  shading  palette  as  used  in  MRI  image 
representation.  The  whole  array  was  then  displayed  as  a 
bitmap,  in  a  manner  similar  to  the  display  of  an  MRI  image.  In 
this  way,  any  algorithms  developed  and  tested  on  the  model 
images,  could  be  directly  applied  to  functional  MRI  images. 

For  the  calculation  of  transformations,  the  familiar  3- 
dimensional  axis  system  was  used  (figure  3)  .  Motion  was 
calculated  as  a  combination  of  rotation  about,  and  translation 
along  the  axes.  Given  3  axes,  a  total  of  6  degrees  of  freedom 
could  be  modelled.  Rotation  and  translation  were  calculated  in 
the  following  forms  : 


P'=P+T 


P/=P'R 


where  P  is  the  point  of  interest  in  vector  form,  and  T  and  R 
are  the  matrix  representations  of  translation  and  rotation 
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respectively.  The  following  matrices  were  used: 
Translation: 

x  Dx 
T  y  =  Dy 
z  Dz 


Rotation  about  z  axis: 

cos  0  sin  0  0 
Rz  ( 0 )  =  -sin  0  cos  0  0 
0  0  1 

Rotation  about  y  axis: 

cos  0  0  -sin  0 
Ry(Q)=  0  10 

sin  0  0  cos  0 

Rotation  about  x  axis: 

10  0 
RX{Q)  =0  cos  0  sin  0 
0  -sin  0  cos  0 

Transformations  to  the  ellipsoid  model  were  calculated  by  the 
application  of  the  above  matrices  to  each  point  of  the  image 
array.  Any  motion  could  be  represented  as  a  combination  of 
translation  and  rotation,  thereby  breaking-down  a  difficult 
transformation  into  a  number  of  simpler  steps.  To  move  an 
object  between  two  locations  in  3  dimensional  space,  the 
following  steps  were  followed: 

1.  shift  of  origin  to  point  of  rotation 

2.  rotate  about  desired  axis 


. 

. 
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3.  translate  along  desired  axis 

4.  move  origin  back  by  original  offset 

By  application  of  the  above  steps  to  every  point  of  the  image 
array,  any  motion  in  space  could  be  performed. 

Intensity  variations  of  the  ellipsoid  were  modeled  to  be 
comparable  to  the  intensity  variations  found  in  a  typical  MR  I 
image.  Based  on  a  palette  of  256  shades  of  gray,  the  model 
uses  a  random  number  generator  to  calculate  intensities  in  the 
realm  of  the  MR I  shading  palette.  A  representation  of  typical 
MRI  intensities  was  thereby  obtained.  Similarly,  an  artificial 
level  of  noise  could  be  added  to  the  model  image.  Though 
levels  of  noise  in  MRI  images  are  often  not  high,  they  are 
still  present,  and  are  a  reflection  of  imperfections  in 
transmission  and  processing  techniques.  Noise  is  seen  as 
pixels  randomly  spread  throughout  the  image.  By  application  of 
a  simple  noise  generation  algorithm63  the  model  allows  the  user 
to  add  a  level  of  noise  equal  to  approximately  6%  of  the  image 
value.  The  various  perturbations  and  functions  allowed  by  the 
model  were  programmed  to  be  interactively  controlled  by  the 
user.  The  computer  display  shows  the  two  images  being  compared 
(MRI  or  model  images)  ,  as  well  as  the  operator  controls 
(figure  6)  .  The  user  may  therefore  manipulate  the  image  at 
will,  and  apply  the  desired  detection  techniques  to  the 
manipulated  image. 
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III.  Perturbed  Echo  Planar  Image  Sets 

After  initial  trial  with  the  ellipsoid  model  described 
above,  the  same  techniques  were  evaluated  with  controlled 
perturbation  of  actual  echo-planar  image  sets.  The  image  sets 
used  were  chosen  from  the  library  of  echo  planar  images 
available  at  the  neuropsychology  laboratory  at  the  West  Haven 
VA  Hospital.  They  were  images  that  had  been  used  in  previous 
functional  echo-planar  studies.  The  images  were  converted  from 
a  256  X  128  pixel  byte-swapped  floating  point  format,  and 
interpolated  to  a  256  X  256  floating  point  image  by  an  MRI 
package  written  for  the  neuropsychology  laboratory  by  Anthony 
Adrignollo.  This  256  X  256  floating  point  array  was  then  read 
by  the  author's  computer  program  for  evaluation  by  the  above 
algorithms. 

Perturbation  of  the  echo-planar  image  in  the  plane  of 
imaging  was  imitated  by  application  of  translation  and 
rotation  equations  in  the  same  manner  as  for  the  ellipsoid 
model.  A  similar  procedure  was  followed  as  with  the  ellipsoid 
model,  with  evaluation  of  the  movement  detection  algorithms 
with  the  echo-planar  image.  Motion  in  the  plane  of  imaging  and 
perpendicular  to  plane  of  imaging  were  again  evaluated. 
Movement  in  the  plane  of  imaging  was  monitored  to  the  nearest 
pixel  on  the  screen.  Given  the  image  resolution  of  256  X  256 
pixels,  with  the  field  of  view  of  40  X  40  cm  produced  by  an 
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echo-planar  image,  each  pixel  correlated  with  1.56mm 
translation  of  the  imaged  object.  Movement  out  of  the  plane  of 
imaging  was  modeled  by  interpolation  between  two  consecutive 
images  along  the  z  axis.  Since  such  images  were  always 
separated  by  7  mm  in  space,  interpolations  were  performed  at 
the  equivalent  of  1  mm  separations. 

TV.  Phantom  Studies 


The  movement-detection  algorithms  were  finally  evaluated 
by  application  to  an  image  set  of  a  phantom  model  of  a  human 
head.  A  3  dimensional  hollow  plastic  model  of  a  human  head  was 
used.  (Figure  7)  To  give  the  interior  of  the  model  an 
appearance  similar  to  a  human  head,  a  rubber  hose  was  placed 
in  the  approximate  location  of  the  sagittal  sinus  and  the 
hollow  model  filled  with  a  gel.  The  phantom  was  mounted  on  a 
plastic  mount  which  would  allow  various  motions  of  the  model. 
By  moving  the  model,  in  combination  with  motion  of  the  gantry 
bed,  motions  in  and  out  of  the  plane  of  imaging  were  imitated. 
Motions  for  translation,  and  rotation  were  performed. 

The  images  were  taken  with  the  standard  echo-planar 
imaging  sequence  that  had  been  used  in  acquisition  of  previous 
echo-planar  studies.  The  images  were  obtained  at  a  resolution 
of  128  X  256  pixels,  with  the  corresponding  data  set 
interpolated  and  processed  in  the  same  manner  as  was  done  with 
the  echo-planar  images  in  the  previous  phase  of  the  study. 
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Figure  6 

Photograph  of  the  program  screen  .  Two  echo- 
planar  images  may  be  compared. 


Figure  7 

Photograph  of  the  phantom  model  used  for  the 
phantom  phase  of  the  study. 
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RESULTS 

All  three  methods  of  movement  detection  were  evaluated 
with  the  model  images,  the  echo-planar  images  and  the  phantom 
images.  For  each  case,  results  for  movement  in  the  plane  of 
imaging  will  be  presented  before  results  for  movement  out  of 
the  plane  of  imaging.  All  figures  may  be  found  at  the  end  of 
the  results  section. 

Center  of  Mass 

Figure  8  shows  the  results  of  the  application  of  the 
center  of  mass  algorithm  to  the  ellipsoid  model,  with  added 
noise  that  had  been  filtered.  Linear  translation  along  the  x 
and  y  axes  could  be  reliably  detected  by  this  approach.  Linear 
regression  analysis  applied  to  the  data  points  showed  a 
correlation  coefficient  R2  of  0.9961,  a  slope  (m)  of  0.9285 
and  an  y0 intercept  of  0.9376.  Rotation  in  the  plane  of  imaging 
of  the  ellipsoid,  however,  did  not  result  in  significant 
change  in  center  of  mass  coordinates,  as  may  be  seen  in  figure 
8b,  with  linear  regression  analysis  reflecting  that  fact. 

The  same  algorithm  was  applied  to  a  set  of  actual  echo- 
planar  images,  as  seen  in  Figure  9.  A  total  of  11  EPI  images 
with  660  total  perturbations  in  the  plane  of  imaging  were 
analyzed  (440  linear  translations  and  220  rotations) .  Again, 
as  with  the  ellipsoid  model,  linear  movement  in  the  plane  of 
imaging  was  found  to  correlate  well  with  changes  in  the  center 
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of  mass  coordinates.  Regression  analysis  showed  a  correlation 
coefficient  of  0.9991  for  linear  translation  in  plane,  with  m 
and  y0  values  of  0.8669  and  -0.0746  respectively.  Center  of 
mass  shift  was  within  2  pixels  of  actual  movement  along  the 
x  and  y  axes  in  93%  of  the  cases  (408/440)  .  The  center  of 
mass  algorithm  did  not  overestimate  the  amount  of  shift,  but 
did  underestimate  the  movement  in  44%  of  the  cases  (193/40). 
Changes  in  center  of  mass  coordinates  in  the  case  of  rotation 
of  the  actual  echo-planar  images,  however,  did  not  correspond 
to  actual  movement  (Figure  9b)  .  Though  shift  occurred  with 
rotation  of  the  images,  it  could  not  be  correlated  with  change 
of  center  of  mass  coordinates. 

Finally,  the  center  of  mass  algorithm  was  applied  to 
movement  in  the  plane  of  imaging  of  the  phantom  (Figure  10) . 
Results  for  linear  translation  in  the  plane  of  imaging  was 
consistent  with  previous  results,  showing  good  correlation, 
with  a  R2  value  of  0.9972,  m  value  of  0.9632  and  a  value  of 
0.1764  for  the  y  intercept.  While  rotation  in  the  plane  of 
imaging  did  result  in  a  small  shift  in  center  of  mass 
coordinates  (center  of  mass  not  0,0),  this  shift  appeared  to 
stay  virtually  constant  within  the  range  of  movement  tested 
(Figure  10b) . 

The  center  of  mass  algorithm  was  then  tested  as  a  method 
of  detection  of  movement  out  of  the  imaging  plane.  In  Figure 
11  one  may  see  the  results  as  found  with  the  ellipsoid  model. 
Change  in  center  of  mass  coordinates  did  occur,  but  was  found 
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not  to  correlate  in  a  meaningful  manner  with  either 
translation  through  plane  along  the  z  axis,  or  rotation 
through  plane  around  the  x  and  y  axes  (Linear  correlation 
calculations  as  given  in  the  figure) . 

Center  of  mass  calculations  performed  on  echo-planar 
images  with  a  total  of  220  perturbations  showed  some  increased 
variability  in  the  center  of  mass  coordinates  than  was  seen  in 
the  ellipsoid  (see  Figure  12) .  However,  as  in  the  case  in  the 
case  of  the  ellipsoid,  no  meaningful  correlation  with  actual 
movement  could  be  obtained. 

This  result  extended  also  to  center  of  mass  shift  in  the 
images  of  the  phantom  model.  Figure  13  shows  the  results  of 
center  of  mass  calculations  as  performed  with  the  phantom. 
Considerable  fluctuations  may  be  seen  in  the  case  of 
translation  through  plane,  with  significant  sharp  changes  in 
the  center  of  mass  coordinates.  As  expected,  linear  regression 
showed  poor  correlation  with  actual  movement.  In  the  case  of 
rotation  in  the  plane  of  imaging,  a  gradual  change  of  the 
center  of  mass  coordinates  was  observed. 

Variance  of  Ratios 

The  variance  of  ratios  algorithm  was  designed 
specifically  for  the  detection  of  movement  in  the  plane  of 
imaging,  and  it  demonstrated  itself  capable  in  this  regard. 
Figure  14  shows  the  results  of  the  algorithm  as  applied  to  an 
ellipsoid  model,  as  well  as  to  perturbed  echo-planar  images. 
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In  both  cases,  perfect  detection  could  be  attained  for  linear 
translation  in  the  plane  of  imaging,  as  well  as  for  rotation 
around  the  z  axis.  This  method  distinguished  itself  from  the 
center  of  mass  approach  by  being  able  to  accurately  detect 
rotation  in  the  plane  of  imaging.  Linear  regression  analysis 
calculated  values  of  R2,  m,  and  y0,  of  1.0,  1.0,  and  0  for  both 
translation  and  rotation  in  the  plane  of  imaging. Phantom 
results  were  similar,  with  very  good  detection  and 
quantification  of  movement  for  both  linear  translation  as  well 
as  rotation  in  the  plane  of  imaging.  The  calculated  R2  value 
for  linear  translation  in  the  imaging  plane  was  0.9991,  the  m 
value  was  0.9877  and  the  y0  value  was  0.0392.  In  the  case  of 
rotation,  very  good  correlation  was  found  over  the  range  of 
movement  tested,  with  linear  correlation  analysis  giving 
values  of  0.9865  for  correlation  coefficient  R2,  1.0545  for 
slope  m,  and  0  for  y-intercept  y0. 

While  this  method  was  designed  only  for  detection  of 
movement  in  the  plane  of  imaging,  an  attempt  was  made  to  get 
an  indication  of  movement  out  of  the  plane  of  imaging  by 
monitoring  of  the  minimum  value  of  the  variance  of  ratios 
converged  to  (v,,^)  .  As  described  in  methods,  the  algorithm 
searches  for  the  perturbation  of  images  with  the  lowest  v 
(variance  of  ratios)  value,  and  designates  that  value  vm.  The 
corresponding  perturbation  is  treated  as  the  detected 
movement.  In  principle,  a  poor  fit  would  result  with  a  higher 
value  of  v^  than  would  a  good  fit.  Values  of  v^  converged  to 
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were  compared  for  images  known  to  be  out-of-plane  with  each 
other,  as  well  as  with  images  known  to  be  in  the  plane  of 


TABLE  1 

Ellipsoid  EPI  Phantom 


v0 

Vi 

v„ 

Vi 

v„ 

Min 

-0.03181 

0.08311 

0 

0.06193 

-0.02183 

0.00837 

Max 

0.05162 

0.14092 

0 

0.11888 

0.03748 

0.38048 

SD 

0.01933 

0.01460 

0 

0.01602 

0.01729 

0.09470 

Mean 

0.0 

0.10918 

0 

0.09247 

0.0 

0.12544 

p  <  1.0 

e  -5 

p 

<  1.0  e  -5 

P  < 

1.0  e  -5 

V;  =  normalized  v,,^  for  two  images  in  the  plane  of  imaging 
v0  =  normalized  vinjn  for  two  images  out  of  the  plane  of 
imaging  to  each  other 


imaging.  This  was  performed  for  the  model  ellipsoid,  the  set 
of  actual  perturbed  echo-planar  images,  as  well  as  for  echo- 
planar  images  obtained  of  the  phantom.  The  results  may  be  seen 
in  table  1.  The  results  were  normalized  by  calculating  the 
average  of  the  values  of  v^  in  the  plane  of  imaging,  and 
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subtracting  that  value  from  all  other  values  of  v,^.  As  may  be 
seen,  results  were  relatively  consistent  over  the  different 
phases  of  the  study,  indicating  the  ability  to  detect  movement 
out  of  the  imaging  plane  by  this  method.  P  values  calculated 
for  the  differences  were  found  to  be  significant  in  all  cases. 

Area 

As  expected,  analysis  by  the  area  method  was  found  to  be 
poorly  predictive  of  movement  in  the  plane  of  imaging.  Figure 
16  shows  the  results  of  the  area  algorithm  applied  to  an 
ellipsoid  with  added  nose  that  had  been  filtered.  Area  changes 
as  a  percentage  of  the  total  area  of  the  original  ellipse  were 
found  to  be  0  for  translation  along  both  the  x  and  y  axes,  and 
ranging  from  -0.6894%  to  2.2487%,  with  an  average  of  0.3286  % 
for  rotation  in  the  plane  of  imaging.  Though  theoretically  one 
would  not  expect  any  change  at  all  in  area  for  rotation  in 
plane,  the  finite  pixel  size  of  the  screen  adds  a  small  noise 
factor  to  area  calculations,  with  a  resultant  variation  in 
area . 

To  get  a  better  understanding  of  the  size  of  this  factor 
so  as  to  be  able  to  differentiate  it  from  variations  in  area 
due  to  movement  through  the  plane  of  imaging,  further  studies 
were  performed  with  echo-planar  images.  For  area  comparisons 
with  perturbed  echo-planar  images,  11  sets  of  images  with  a 
total  of  220  in-plane  rotations  were  studied.  Small  changes 
were  noted  with  rotations  in  plane,  ranging  from  0  to  3.6239%, 
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with  the  average  value  being  given  by  0.47706%  (figure  18a). 
Echo-planar  images  of  the  phantom  rotated  in  the  plane  of 
imaging  also  showed  even  smaller  variations  in  area  size, 
ranging  from  0  to  0.5797%,  with  an  average  of  0.3087%  (figure 
16b)  . 

The  area  method  did  prove  to  be  effective  in  detection  of 
movement  out  of  imaging  plane.  Figure  17  shows  the  percentage 
area  change  with  movement  perpendicular  to  the  plane  of 
imaging  for  ellipsoids  with  added  nose  that  had  been  filtered. 
The  rate  of  change  was  seen  to  be  non-linear,  with  further 
dependence  on  the  thickness  of  the  object  being  moved  (c  term 
in  ellipsoid  eguation) .  As  expected,  a  more  rapid  decrease  of 
area  is  found  for  ellipsoids  of  lesser  thickness.  Variations 
in  area  for  rotation  of  the  ellipsoid  model  through  the  plane 
of  imaging  were  seen  to  be  consistently  less  than  <  2.5%,  with 
no  particular  pattern  seen  with  the  movement.  The  variations 
were  no  larger  than  the  variations  in  area  seen  with  rotations 
in  the  plane  of  imaging. 

Application  of  the  area  algorithm  to  perturbed  echo- 
planar  images  confirmed  the  pattern  (Figure  18) .  Eleven  echo- 
planar  image  sets  with  a  total  of  77  out-of-plane 
interpolations  were  compared.  As  may  be  seen  in  the  figure  and 
table  2,  the  area  change  was  variable,  ranging  from  a  minimum 
of  0.04%  to  17.224  %  of  the  total  area. 

Area  analysis  of  the  echo-planar  images  of  the  phantom 
showed  large  changes  in  area  (figure  19)  .  Knowing  that  the 
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changes  would  depend  on  the  segment  of  the  head  being  imaged. 


TABLE  2 


Distance  Hin 


Max  Mean  St. Dev 

-2.09085  -0.09590  0.20268 
-4.59241  -2.09573  1.52406 
-7.66646  -3.51221  2.64182 
-12.0473  -5.84649  4.11626 
-14.2004  -6.97743  4.67183 
-16.0921  -8.41335  5.19402 
-17.2246  -9.26353  5.42246 


1  -0.04128 

2  -0.09841 

3  -0.07567 

4  -0.18575 

5  -0.20638 

6  -0.37149 

7  -0.54348 


two  different  regions  of  the  phantom  were  imaged  for  detection 
of  translation  through  the  plane  of  imaging.  Area  change  was 
monitored  towards  the  base  of  the  head,  and  towards  the  apex, 
where  greater  changes  were  expected.  Changes  as  large  as  80% 
were  observed  for  a  movement  of  12  mm  towards  the  apex  (mean 
change  25.5%).  For  rotation  through  the  plane  of  imaging, 
variations  in  area  were  found  to  be  smaller,  but  still 
significant,  ranging  from  3.3%  to  19.2%,  with  a  mean  of  10.7%. 
Larger  changes  in  area  were  therefore  observed  with  rotation 
through  the  plane  of  imaging  of  the  phantom  head  than  with  the 
ellipsoid  model,  indicating  the  ability  to  be  able  to  detect 
such  rotating  movement. 
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actual  movement  (pixels) 


Fig  8a 


Fig  8b 


Figure  8 

Figure  of  the  center  of  mass  algorithm  applied 
to  an  ellipsoid  with  a=60,b=80,c=30,  with  added 
noise  that  had  been  filtered,  imaged  through 
the  plane  z=0.  Figure  8a  shows  the  center  of 
mass  change  as  a  function  of  linear  translation 
along  both  x  and  y  axes.  Figure  8b  shows  the 
change  in  center  of  mass  coordinates  for 
rotation  in  the  plane  of  imaging. 
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Fig  9a 


— .  10  1 

« 

03 

.a  5  - 

y 

1  t  |  | 

+  *+  +  +  4-  +  +■ 

3  X  +  +  -5  . 

o  ‘r 

°  -10  J 

-10  -5  C 

actual  rotatic 

■  f 

X  u  8 

)  5  10 

>n  (degrees) 

Fig  9b 


Figure  9 

Center  of  mass  algorithm  applied  to  11  echo- 
planar  images  with  a  total  of  660 
perturbations  in  the  plane  of  imaging.  Figure 
9a  shows  change  of  center  of  mass  coordinates 
with  linear  transformation.  Figure  9b  shows 
center  of  mass  change  with  rotation  in  the 
plane  of  imaging. 
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Figure  10 

Center  of  mass  algorithm  applied  to  echo-planar 
images  obtained  of  the  phantom  model  moved  by 
controlled  amounts  in  the  plane  of  imaging.  Figure 
10a  shows  the  results  with  linear  movement  in  the 
plane  of  imaging.  Figure  10b  shows  the  results  of 
rotation  in  the  plane  of  imaging. 
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Fig  11a 


Fig  lib 


Figure  11 

Results  of  the  center  of  mass  algorithm  applied 
to  an  ellipsoid  with  a=60, b=80,c=30  imaged  in  the 
plane  z=0.  Figure  11a  shows  the  center  of  mass 
coordinates  for  successive  cross-sectional  images 
of  the  ellipsoid  as  it  moves  along  the  z-axis. 
Figure  10b  shows  the  center  of  mass  coordinates 
for  images  of  the  ellipsoid  as  it  rotates  through 
the  plane  of  imaging. 

Fig.  11a  Fig.  lib 
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Fig  12a 


Fig  12b 


Figure  12 

Center  of  mass  algorithm  applied  to  a  set  of 
11  echo-planar  images  with  a  total  of  77  out 
of  plane  interpolations.  Linear  interpolation 
along  the  z  axis  were  performed.  Figure  12a 
shows  the  change  in  center  of  mass  x 
coordinate  with  movement  along  the  z  axis, 
while  Figure  12b  shows  the  change  in  y 


coordinate. 
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Fig  13a 


Fig  13b 


Figure  13 

Results  of  the  center  of  mass  algorithm  applied 
to  echo-planar  images  of  the  phantom  model. 
Figure  13a  shows  center  of  mass  change  for 
translation  along  the  z  axis.  Figure  13b  shows 
change  in  center  of  mass  coordinates  for 
rotation  of  the  phantom  through  the  plane  of 
imaging. 

Figure  13a  Figure  13b 
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Fig  14a 


Fig  14b 


Figure  14 

Variance  of  ratios  algorithm  as  applied  to  an  ellipsoid 
with  a=60,b=80,c=30  imaged  through  the  plane  z=0,  as 
well  as  a  set  of  echo  planar  images  with  a  total  of  660 
perturbations.  Figure  14a  shows  the  detected  movement 
as  a  function  of  linear  translation  in  the  plane  of 
imaging,  while  figure  14b  shows  detected  movement  as  a 
function  of  rotation. 
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Fig  15a 


Figure  15 

Variance  of  ratios  algorithm  as  applied  to  echo- 
planar  images  of  the  phantom  model  for  movement  in 
the  plane  of  imaging.  Figure  15a  shows  the  detected 
linear  translation  as  a  function  of  actual  movement 
along  the  x  and  y  axes.  Figure  15b  shows  the 
detected  rotation  as  a  function  of  actual  rotation 
in  the  plane  of  imaging. 
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Fig  16a 


Fig  16b 


Figure  16 

Area  change  algorithm  as  applied  to  the  model 
ellipsoid  and  echo-planar  images  of  the  phantom 
model.  (Results  for  set  of  echo-planar  images  may  be 
seen  in  figure  18a) .  Figure  16a  shows  area  change 
with  linear  translation  and  rotation  in  the  plane  of 
imaging.  Figure  16b  shows  area  change  with  rotation 
of  the  phantom  model  in  the  plane  of  imaging. 
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Figure  17 

Area  algorithm  applied  to  the  ellipsoids  of  various 
sizes  with  added  noise  that  had  been  filtered.  Figure 
17a  shows  the  changes  in  cross-sectional  area  of 
ellipsoids  with  various  c  values  (60-30)  as  a 
function  of  movement  along  the  z  axis.  Figure  17b 
shows  the  changes  in  cross-sectional  area  of 
ellipsoids  with  various  c  values  (80-50)  as  a 
function  of  rotation  through  the  plane  of  imaging. 
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Fig  18a 


Fig  18b 


Figure  18 

Results  of  the  area  algorithm  as  applied  to  a  set 
of  11  echo-planar  images.  Figure  18a  shows  the 
area  change  for  rotation  in  the  plane  of  imaging. 
Figure  19b  shows  the  results  of  area  change  for 
images  interpolated  to  lie  at  different  points 
along  the  z  axis. 


Page  58 


Fig  19a 


Fig  19b 


Figure  19 

Results  of  area  change  algorithm  as  applied 
to  echo-planar  images  of  the  phantom  model 
for  transformations  through  the  plane  of 
imaging.  Figure  19a  shows  movement  along  the 
z  axis,  while  Figure  19b  shows  rotation 
through  the  plane  of  imaging. 
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DISCUSSION 

Three  different  methods  were  developed  and  tested  in  the 
attempt  to  detect  patient  movement  in  functional  MRI  imaging. 
All  were  found  to  be  useful  in  detecting  some  aspect  of 
motion.  The  methods  will  be  discussed  separately,  with  a 
combined  analysis  to  follow.  In  particular,  each  methods  will 
be  discussed  in  relation  to  its  ability  to  detect  and 
quantitate  the  movement  of  the  imaged  subject. 

Center  of  Mass 

The  center  of  mass  technique  proved  itself  reliable  in 
the  detection  of  movement  in  general.  Based  on  a  physical 
principle,  the  center  of  mass  is  relatively  straight-forward 
to  compute,  and  applicable  to  any  image  set.  The  computation 
itself  is  rapid,  and  does  not  necessarily  require  any  pre¬ 
processing  of  the  image.  As  may  be  seen  from  the  study 
results,  the  center  of  mass  of  an  image  changes  with  virtually 
all  types  of  motion.  In  all  phases  of  the  study,  and 
particularly  with  analysis  of  actual  echo-planar  images, 
changes  in  center  of  mass  coordinates  were  noted  with  movement 
in  plane,  and  out  of  plane.  This,  of  course  is  to  be  expected 
if  one  considers  the  nature  of  the  calculation.  Being  a 
function  of  the  intensity  and  position  of  2562  points,  even 
small  changes  in  position  or  intensity  as  a  result  of  motion 
would  register  as  a  shift  in  the  center  of  mass.  However,  not 
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in  all  cases  was  the  shift  mathematically  predictable  or 
reproducible.  Therefore,  while  this  method  could  detect  most 
movements,  it  could  not  quantitate  all  movements. 

The  predictability  and  reproducibility  of  the  results  is 
key  to  the  application  of  a  technique  for  quantification  of 
motion.  For  example,  in  the  case  of  rotation  in  the  imaging 
plane  of  typical  EPI  images  (Fig  9b) ,  changes  in  both  the  x 
and  y  coordinates  of  the  center  of  mass  could  be  seen  - 
especially  at  greater  rotations.  Similar  changes,  on  a  smaller 
scale,  were  seen  with  rotation  of  the  modelled  ellipsoid  and 
the  phantom  (Fig.  8b,  10b) .  In  all  three  cases,  however,  the 
shifts  did  not  correlate  in  a  predictable  fashion  with  actual 
motion.  Similar  results  were  observed  with  motion  out  of 
plane.  Any  shifts  in  center  of  mass  that  occurred  could  not  be 
reliably  correlated  with  the  actual  movement. 

This  may  be  understood  for  both  rotation  in  plane  and 
movement  out  of  plane.  For  rotation  in  plane,  a  change  in  the 
center  of  mass  coordinates  would  depend  on  the  distance  of  the 
center  of  mass  of  the  object  from  the  point  of  rotation.  If 
the  point  of  rotation  coincides  with  the  center  of  mass  point, 
the  center  of  mass  will  not  shift  with  rotation.  The  center  of 
mass  made  only  a  minor  shift  for  rotation  in  the  plane  of 
imaging  of  the  phantom  model,  for  example  (Figure  10b) .  The 
shift  will  depend  on  the  sine  and  cosine  of  the  angle  of 
rotation  and  the  distance  of  the  point  of  rotation  form  the 
center  of  mass.  With  greater  distance  of  the  center  of  mass 
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from  the  point  of  rotation,  larger  shifts  of  center  of  mass 
coordinates  will  result.  In  the  most  ideal  situation,  an  a- 
priori  knowledge  of  the  center  of  mass  coordinates  relative  to 
the  point  of  rotation  is  needed  for  interpretation  of  center 
of  mass  shift  with  rotation.  As  this  is  not  available 
generally,  correlation  between  center  of  mass  and  rotation  in 
plane  is  not  predictable. 

In  the  case  of  movement  out  of  the  plane  of  imaging,  a 
new  image  is  created,  with  different  anatomical  information, 
and  therefore  different  pixel  intensities.  A  center  of  mass 
shift  will  therefore  occur  as  a  result.  As  the  pixel 
intensities  and  positions  of  the  new  image  are  not  generally 
related  to  those  of  the  old  image  in  a  predictable  fashion, 
the  change  in  center  of  mass  coordinates  is  also  not 
predictable.  As  in  the  case  with  rotation  in  the  plane  of 
imaging,  though  a  shift  in  center  of  mass  occurs,  it  is  not 
possible  to  correlate  this  shift  with  the  degree  of  motion 
through  plane. 

Linear  translation  in  the  imaging  plane  is  the  movement 
that  is  most  reliably  detected  and  quantitated  by  the  center 
of  mass  method.  This  was  found  to  be  the  case  with  the 
ellipsoid  model,  the  echo-planar  images  and  the  phantom  study. 
In  all  situations  a  linear  correlation  between  the  amount  of 
motion,  and  the  center  of  mass  shift  could  be  observed.  This 
method,  therefore  could  provide  a  reliable  estimate  of 
movement  in  the  plane  of  imaging.  To  be  able  to  do  so, 
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however,  the  information  provided  by  this  method  must  be 
placed  into  context,  by  other  data.  A  change  in  center  of  mass 
coordinates  may  be  due  to  a  number  of  motions,  as  discussed 
above.  However  only  if  the  change  is  due  to  a  linear 
transformation  in  plane  does  the  shift  have  any  quantitative 
value.  In  other  words,  further  information  is  needed  regarding 
the  motion  to  be  able  to  properly  interpret  the  results  of  the 
center  of  mass  shift.  This  may,  in  principle,  be  obtained  from 
the  other  methods  of  movement  detection. 

It  should  be  noted,  however,  that  even  in  the  case  of 
linear  movement  in  the  imaging  plane,  only  an  estimate  of 
movement  is  obtained.  The  method  was  not  able  to  detect  linear 
movement  exactly,  though  it  did  estimate  it  to  within  2  pixels 
in  over  90%  of  cases  with  analysis  of  echo-planar  images. 
Furthermore,  the  method  had  a  tendency  to  underestimate  the 
actual  linear  movement.  As  such,  while  this  technique  will 
give  a  reasonable  estimate  of  linear  movement  in  the  imaging 
plane,  in  approximately  half  the  cases  the  technique  will 
underestimate  the  actual  size  of  movement. 

Overall,  this  is  a  quick  and  simple  method  of  movement 
detection.  It's  ease  of  application  and  speed  of  calculation 
make  it  a  good  initial  test  of  movement.  The  method  has  the 
potential  to  detect  all  forms  of  movement.  The  method  cannot, 
however  determine  what  type  of  movement  resulted  in  the  center 
of  mass  shift.  If  , however,  the  movement  is  known  to  be  in 
the  plane  of  imaging,  the  center  of  mass  technique  will  give 
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a  good  reliable  estimate  of  linear  translation.  It  may  thus 
detect  all  movement,  and  approximately  guantify  linear 
movement  in  the  plane  of  imaging. 

Variance  of  Ratios 

The  technique  of  variance  of  ratios  as  implemented  in 
this  project,  is  based  on  the  technique  described  by  Woods 
et.al.  In  general,  it  was  found  to  be  very  reliable  in 
detection  of  motion  in  the  imaging  plane.  Because  it  functions 
by  doing  a  pixel  by  pixel  comparison  of  the  images  being 
analyzed,  it  also  results  in  greater  accuracy.  For  motion  of 
the  ellipsoid  model,  as  well  as  the  perturbed  echo-planar 
images,  the  variance  of  ratios  technique  was  able  to  detect 
and  quantitate  both  translation  as  well  as  rotation  in  plane 
(Fig  14,15). 

As  described  in  the  results  reaction,  an  attempt  was  also 
made  to  use  the  minimum  value  of  the  variance  of  ratios 
converged  upon  (v^)  to  give  an  indication  of  motion  out  of 
plane.  For  perfectly  aligned  images,  the  value  of  v^  should 
be  0,  as  was  the  case  with  the  perturbed  EPI  images.  Even  if 
not  0,  however,  the  transformation  with  the  lowest  v  is 
considered  by  the  algorithm  to  be  the  calculated  movement  in 
plane.  Images  out  of  plane  with  each  other,  could  not  achieve 
a  vmin  value  of  0,  and  the  value  of  v^  converged  upon  for  these 
images  should  be  higher  than  that  of  images  in  plane.  As  the 
results  demonstrate,  (table  1)  this  approach  did  turn-out  to 
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be  effective  in  detection  of  movement  outside  of  the  plane  of 
imaging.  While  this  approach  does  not  quantitate  the  motion 
out  of  plane,  it  may  be  used  as  an  important  indicator  of  such 
motion. 

The  variance  of  ratios  algorithm,  though  a  powerful 
method  for  quantification  of  movement  in  the  plane  of  imaging, 
is  very  calculation  intensive,  and  therefore  slow.  A  value  of 
v  must  be  calculated  for  each  possible  transformation.  In 
search  of  a  best  fit,  the  algorithm  calculates  all  possible 
transformations  of  the  image  in  a  certain  area  of  the  screen, 
recalculating  the  variance  of  ratios  each  time.  To  find  the 
best-fit  in  a  10  by  10  pixel  area  of  the  screen,  a  total  of 
103  transformations  are  needed,  for  each  of  which  a  v  value 
must  be  calculated.  This  corresponds  to  a  total  of 
approximately  5.1  X  1010  calculations  for  this  area  of  the 
screen.  As  a  comparison,  to  determine  the  centers  of  mass  of 
two  images,  approximately  5.2  X  105 calculations  are  needed. 
Clearly,  while  the  variance  of  ratios  algorithm  is  very 
effective  in  detection  and  quantification  of  motion  in  plane, 
it  is  also  very  time  consuming. 

Overall  this  is  a  powerful  technique  for  detection  as 
well  as  quantification  of  motion  in  the  plane  of  imaging.  It 
may  reliably  quantify  linear  motion  in  plane,  as  well  as 
rotation  in  plane.  However,  the  variance  of  ratios  technique 
is  extremely  calculation  intensive,  and  therefore  limited  by 
its  slow  speed  of  operation. 
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Area 

The  area  algorithm  is  effective  in  detection  of  motion 
outside  of  the  plane  of  imaging.  As  described  in  the  methods 
section,  it  uses  the  image  processing  techniques  of  low  pass 
filters,  Laplacian  operators  and  thresholding  to  outline  the 
edges  of  an  image  and  calculate  the  area  of  the  outlined  form. 
As  expected,  this  method  was  not  effective  in  detection  of 
motion  in  the  plane  of  imaging.  In  cases  of  motion  in  the 
plane  of  imaging,  area  changes  were  either  0  (linear 
translation)  or  relatively  small  (rotation  in  plane) .  Because 
the  area  approach  is  based  on  an  outline  and  fill  technique, 
its  accuracy  is  limited  by  the  resolution  of  the  screen.  As 
such,  small  variations  of  area,  corresponding  to  the  finite 
size  of  the  screen  pixel  would  be  expected  even  for  motions  in 
plane.  In  all  cases  this  change  was  found  to  be  a  fraction  of 
one  percent  on  average. 

Dependent  on  changes  in  area  for  detection  of  movement 
out  of  the  plane  of  imaging,  this  technique  relies  upon  the 
asymmetry  of  the  imaged  object.  Therefore,  the  amount  of  area 
change  for  motion  perpendicular  to  the  plane  of  imaging  was 
dependent  on  the  segment  of  the  object  being  imaged.  In  the 
case  of  the  perfect  ellipsoid  model,  the  rate  of  change  was 
greater  for  ellipsoids  of  lesser  thickness,  and  smaller  for 
ellipsoids  of  greater  thickness.  This  of  course  is  intuitive, 
as  different  cross-sections  of  an  object  with  large  variations 
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in  shape  would  show  significant  variations  in  area.  Unlike  the 
ellipsoid,  whose  geometrical  shape  and  behavior  is  well 
defined,  the  changes  in  areas  with  real  images  are  somewhat 
less  predictable.  A  general  pattern  of  change  may  be  observed, 
without  a  definite  mathematical  relationship.  The  changes 
would  be  dependent  on  the  segment  of  the  head  being  imaged. 
For  example,  motion  perpendicular  to  the  plane  of  imaging  is 
easier  to  detect  for  images  in  the  most  anterior  aspects  of 
the  frontal  lobes,  due  to  the  more  rapid  tapering  of  the  brain 
surface  in  that  region  .  On  the  other  hand,  images  taken 
further  posteriorly,  -  coronal  cuts  in  the  region  of  the  brain 
stem  for  example  -  would  show  little  variation  in  area  with 
motion  through  the  plane  of  imaging.  In  the  case  of 
interpolated  echo-planar  images  (Figure  18b) ,  the  largest 
variation  in  area  was  seen  in  cross-sections  of  images  of  the 
anterior  aspect  of  the  frontal  lobe. 

In  the  case  of  rotation  through  the  plane  of  imaging, 
area  changes  were  found  to  be  smaller  than  for  translation  out 
of  plane.  As  in  the  case  of  translation  through  the  plane  of 
imaging,  the  rate  of  change  of  area  depends  on  the  segment  of 
the  imaged  object.  While  area  changes  for  rotation  of  the 
ellipsoid  model  were  small  (Figure  17b) ,  larger  changes  were 
seen  for  rotation  through  the  plane  of  imaging  of  the  phantom. 
With  the  phantom  model,  area  changes  as  large  as  19%  for  a  3 
degree  rotation  were  seen.  Variations  on  this  scale  are  larger 
than  those  seen  for  rotation  in  the  plane  of  imaging,  and  thus 
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could  be  of  use  in  detection  of  movement  out  of  the  plane. 

The  method  of  area  calculation  by  outlining  and  filling 
by  image  processing  techniques  depends  on  a  relatively  high 
degree  of  human  interaction.  The  use  of  low  pass  filters  in 
combination  with  thresholding  and  Laplacians  filters  is  not 
always  successful  in  effectively  outlining  the  area  of 
interest.  The  thresholding  level  necessary  for  complete 
outlining  of  an  image  requires  experimentation  by  the 
operator.  In  the  case  of  this  study,  for  example,  it  was  found 
that  a  low  pass  filter,  with  thresholding  at  90%  of  the 
maximum  value  of  intensity  followed  by  application  of  a 
Laplacian  filter,  worked  well  for  outlining  of  most  model 
ellipsoid  images.  On  the  other  hand,  echo-planar  images  were 
better  outlined  by  application  of  a  low  pass  filter  alone, 
with  subsequent  thresholding  at  50%  of  the  maximum  intensity 
value.  However,  even  with  this  approach,  the  outline 
calculated  for  10  of  220  echo-planar  images  (4.5%)  was  not 
complete  enough  for  application  of  the  fill  procedure.  The 
user  must  therefore  observe  the  procedure  to  determine  whether 
the  results  are  reliable.  This  high  level  of  human  interaction 
necessary  for  the  proper  application  of  the  area  algorithm 
limits  the  capacity  to  automate  this  procedure. 

In  general,  the  area  method  is  an  easy  to  understand, 
relatively  quick  method  to  apply  that  is  effective  for 
detection  of  movement  out  of  the  plane  of  imaging.  Though  pure 
rotation  through  the  plane  of  imaging  generally  results  in 
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smaller  area  changes  than  translation  through  the  plane  of 
imaging,  it  may  still  be  detected.  This  method,  however  is 
limited  by  the  relatively  high  level  of  user  interaction 
needed  to  make  it  run,  as  well  as  its  inability  to  quantitate 
movement  through  the  plane  of  imaging. 

Combination  of  Techniques 

Each  of  the  above  analyzed  techniques  has  strengths  for 
particular  applications,  as  well  as  limitations  of  its 
capabilities.  A  combination  of  these  techniques  would  be 
helpful  in  detection  of  a  wider  spectrum  of  motion.  Table  3 
shows  the  characteristics  of  the  3  different  techniques.  A 


TABLE  3 


DETECTION 


QUANTIFICATION 


C  Of  M 


all  movement 


linear 

translation  in 
the  imaging  plane 


Var  Ratios 


all  movement 


movement  in  the 
imaging  plane 


Area 


movement  out 
of  the  imaging 
plane 


none 


while  the  variance  of  ratios  technique  can,  in 
principle,  detect  movement  in  and  out  of  the  plane  of 
imaging,  its  slow  speed  limits  its  potential  as  a 
screening  test  for  movement 


. 


Page  69 

proposed  flow-sheet  showing  the  possible  combination  of  these 
techniques  may  be  seen  in  figure  20. 

In  determining  the  order  of  the  steps  of  the  flow-sheet, 
the  strengths  and  weaknesses  of  each  technique  were 
considered.  Based  on  the  abilities  of  the  algorithm,  they  are 
each  assigned  a  certain  task.  In  effect,  therefore,  three 
issues  are  dealt  with:  a)  detection  of  motion;  b) 
classification  of  the  motion  as  being  either  in  or  out  of  the 
plane  of  imaging;  c)  if  in  the  plane  of  imaging,  determination 
the  direction  of  the  motion.  These  three  issues  are  dealt  with 
by  the  center  of  mass  method,  the  area  method,  and  the 
variance  of  ratios  method  respectively. 

The  center  of  mass  technique,  due  to  its  speed  and 
simplicity,  as  well  as  its  ability  to  detect  motion  in 
general,  is  the  initial  test  applied.  If  it  shows  no  variation 
in  center  of  mass  coordinates,  the  probability  of  motion 
occurring  may  be  assumed  to  be  very  small.  By  application  of 
this  technique,  an  initial  estimate  of  motion  may  also  be 
obtained. 

As  the  center  of  mass  technique  does  not  by  itself  give 
any  indications  as  to  the  nature  of  the  motion,  the  area 
algorithm  is  applied  next.  This  puts  the  information  gained  by 
the  center  of  mass  approach  into  context,  by  answering  whether 
the  detected  motion  was  most  likely  in  the  plane  of  imaging  or 
not.  A  significant  change  in  area  indicates  a  motion  through 
the  plane  of  imaging,  while  no  change,  or  minor  change. 
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implies  motion  in  the  plane  of  imaging.  The  cut-off  point  set 
for  detection  of  movement  out  of  the  plane  of  imaging  would  be 
based  on  the  percentage  area  change  observed  in  calculating 
areas  of  echo-planar  images  rotated  in  the  imaging  plane. 

Finally,  if  the  area  algorithm  indicates  a  motion  in  the 
plane  of  imaging,  the  variance  of  ratios  algorithm  may  be 
applied  to  determine  the  most  likely  motion  in  the  plane.  The 
initial  estimates  of  motion  may  be  obtained  from  the  results 
of  the  center  of  mass  algorithm,  with  the  changes  in  center  of 
mass  coordinates  being  input  as  the  starting  linear 
translations  for  the  variance  of  ratios  method.  The  area  of 
the  screen  examined  in  this  situation  could  be  limited  to 
within  a  few  pixels  of  the  initial  estimate,  within  the  degree 
of  error  provided  by  the  center  of  mass  technigue.  The  value 
of  converged  upon  by  this  method  would  be  used  as  a 
confirming  factor  for  movement  in  the  imaging  plane.  Adjusted 
values  of  v^  (as  described  in  results)  below  a  set  cut-off 
would  be  indication  of  a  good  fit,  while  values  above  the  cut¬ 
off  would  indicate  a  poor  fit,  and  interpreted  as  an 
indication  of  movement  out  of  the  imaging  plane. 


Further  Issues 

Each  technique  has  limitation  and  weaknesses  as  already 
discussed.  Similarly,  the  combination  of  all  three  also  has 
significant  limitations  in  movement  detection. 


* 
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In  developing  and  testing  these  techniques,  movement  was 
analyzed  from  the  perspective  of  movement  in  the  plane  of 
imaging  and  out  of  the  plane  of  imaging.  This  approach  was 
chosen  to  simplify  the  problem  into  two  more  easily  manageable 
problems.  The  techniques  tested  are  a  reflection  of  this 
approach.  While  there  is  cross-over  in  the  applicability  of 
each  technique  to  the  two  different  aspects  of  the  problem, 
each  technique  is  better  at  detection  of  one  of  the  two  forms 
of  movement.  No  method  can  reliably  determine  by  itself  a 
combination  of  the  two  motions,  and  the  combination  of  the 
methods  does  not  do  much  better  at  detecting  such  movement. 
Movement  both  through  the  plane  of  imaging  and  in  the  plane  of 
imaging  would  be  classified  as  movement  through  the  plane  by 
the  application  of  the  methods  as  described  by  the  flow-sheet. 

The  heart  of  this  weakness  lies  in  the  inability  to 
quantify  movement  out  of  the  plane  of  imaging.  This  subset  of 
the  problem  is  considerably  more  difficult  than  that  of 
movement  in  the  plane.  In  particular,  the  difficulty  lies  in 
the  lack  of  information  of  what  may  lie  out  of  the  imaging 
plane.  To  be  able  to  make  reliable  predictions  of  movement,  a 
variable  is  monitored  whose  behavior  with  the  change  of 
interest  is  known.  For  movement  in  the  plane  of  imaging,  much 
information  is  available  about  the  contents  of  the  image,  and 
any  perturbations  it  may  undergo.  The  variance  of  ratios 
technique,  for  example,  uses  data  from  every  pixel  of  the 
image  in  its  analysis.  In  applying  a  method  like  the  area 
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method,  an  attempt  is  made  to  use  a  variable  with  some 
predictive  value  over  different  imaging  planes.  However,  while 
it  gives  an  indication  of  movement  out  of  plane,  it  has  no 
ability  to  predict  the  exact  direction  or  combination  of 
movements.  With  movement  in  the  plane  on  the  other  hand,  where 
all  the  image  information  is  available  in  principle,  a  better 
estimate  of  the  direction  of  movement  is  possible  (both  center 
of  mass  and  variance  of  ratios  give  estimates) . 

One  solution  to  this  problem  would  be  the  acquisition  of 
multiple  echo-planar  images  at  small  intervals  throughout  the 
imaging  volume.  By  acquiring  such  an  image  set,  algorithms 
such  as  center  of  mass  and  variance  of  ratios  could  be  applied 
to  the  whole  3  dimensional  set.  Currently  at  most  4  images  at 
7mm  out-of-plane  intervals  are  available  for  comparison.  In 
future,  however,  as  echo-planar  sequences  allow  for  higher 
numbers  of  images  for  comparison  (1mm  cuts  through  the  whole 
volume) ,  better  movement  detection  for  movement  out  of  the 
imaging  plane  may  become  available. 

As  it  stands,  however,  the  limitations  posed  by  the  above 
algorithms  are  quite  serious.  While  movement  may  be  detected 
in  general,  it  can  only  be  quantified  approximately  in  the 
imaging  plane,  and  not  at  all  outside  the  imaging  plane.  In  a 
realistic  scenario,  where  movement  out  of  the  imaging  plane  is 
very  likely,  this  would  amount  to  detection  of  movement,  with 
no  ability  to  be  able  to  determine  its  direction  or 
significance.  Serious  limitations  are  therefore  placed  on  the 
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applicability  of  the  algorithms. 

Given  the  current  limitations,  the  question  arises  as  to 
what  may  be  done  with  the  information  regarding  possible 
movement.  Two  possibilities  may  be  considered.  An  attempt  may 
be  made  to  correct  the  data  set  for  movement,  or  the  data  set 
may  be  deemed  invalid  for  further  analysis  and  discarded. 

If  movement  is  determined  to  be  in  the  plane  of  imaging, 
the  images  could  relatively  easily  be  corrected  by  perturbing 
them  by  application  of  translation  and  rotation  matrices,  as 
discussed  in  theory  and  methods.  Correction  in  this  case  is 
possible  as  all  the  data  necessary  for  proper  reconstruction 
are  available.  Furthermore,  the  correction  would  most  likely 
be  accurate  if  the  variance  of  ratios  algorithm  is  used  to 
determine  the  movement. 

Movement  out  of  imaging  plane,  however,  would  pose  a 
problem  regarding  correction.  Only  an  indication  of  motion 
through  the  plane  of  imaging  may  be  obtained  by  application  of 
the  tested  algorithms.  Since  a  direction  and  quantification  of 
movement  is  not  available,  correction  for  such  movement  would 
not  be  possible.  More  approximate  methods  of  correction  could 
be  considered.  For  example,  one  could  attempt  a  rough 
correction  by  simple  mapping  of  the  image  with  movement  unto 
the  baseline  image.  This,  in  effect  would  amount  to  assuming 
that  the  motion  through  the  plane  of  imaging  was  due  to  a 
translation  through  the  plane  of  imaging,  with  an 
interpolation  being  performed  as  a  correction  factor.  This 
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approach  would  amount  to  an  approximation  at  best,  and  could 
possibly  lead  to  artifacts  itself. 

However,  even  if  no  correction  is  possible,  it  may  be 
argued  that  the  pure  detection  of  movement  is  very  valuable  in 
itself.  At  the  very  least  it  would  alert  the  experimenter  of 
problems  with  the  data  set,  and  allow  him  or  her  to  dispose  of 
the  data  as  desired.  With  the  application  of  the  above 
methods,  this  scenarios  would  arise  for  movement  out  of  the 
plane  of  imaging.  The  decision  to  dispose  of  the  data  set 
would  naturally  depend  on  numerous  factors.  One  could  be  the 
approximate  size  of  the  movement.  While  the  area  detection 
algorithm  does  not  quantitate  the  movement,  it  may  give  an 
approximation  of  the  severity  of  the  movement  in  the  form  of 
the  size  of  area  change.  A  small  change  may  be  deemed 
acceptable,  and  the  data  processed  with  that  fact  in  mind.  A 
further  consideration  could  be  the  difficulty  of  repeating  the 
study.  Repetition  of  the  study  would  be  easier  in  the  case  of 
an  experimental  protocol  than  with  a  patient  having  the  study 
performed  for  diagnostic  purposes. 

Given  the  limitations  of  the  tested  algorithms,  other 
methods  of  dealing  with  the  problem  of  patient  movement  must 
be  considered.  The  tested  algorithms  limited  themselves  to  the 
use  of  raw  echo-planar  images.  The  application  of  an  external 
fiduciary  system,  while  introducing  greater  complexity  into 
the  imaging  process  and  losing  the  capacity  of  retrospective 
application  could  be  useful  in  better  quantification  of 
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movement  through  the  imaging  plane.  More  complete  3 
dimensional  echo-planar  data  sets,  as  they  become  available, 
would  allow  for  the  application  of  the  variance  of  ratios 
algorithm  in  3  dimension,  as  discussed  above.  Finally,  better 
restraining  system  may  be  necessary  as  the  technigue  of 
functional  echo-planar  imaging  becomes  more  widely  clinically 
used. 

As  discussed  in  the  introduction,  the  technique  of 
functional  MRI  imaging  is  a  new  approach  to  functional 
imaging.  Due  to  different  acquisition  methodology  and  much 
slower  acquisition  time,  the  issue  of  movement  detection 
between  images  has  not  been  of  the  same  concern  with  older 
methods  of  MRI  imaging.  This  is  a  problem  that  is  particular 
to  functional  echo-planar  imaging  itself,  and  as  of  the  time 
of  this  writing,  no  previous  studies  dealing  with  movement 
detection  in  functional  MRI  imaging  by  post-processing 
techniques  have  been  published  in  the  literature.  A  study  of 
theoretical  aspects  of  movement  in  echo-planar  imaging  by 
Duerk66  examined  the  ability  to  limit  movement  artifact  in 
echo-planar  imaging  by  modification  of  the  echo-planar 
acquisition  technique  itself.  He  found  that  he  could  decrease 
the  amount  of  signal  loss,  ghosting  and  blurring  in  images  by 
development  of  a  new  theoretical  EPI  acquisition  sequence 
which  refocused  motion  derivatives  at  echoes  and  changed  the 
acquisition  through  moment-nulling  waveforms.  No  studies  have 
addressed  the  issue  of  movement  detection  in  images  that  have 
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already  been  acquired,  however.  This  study  is  therefore  unique 
in  taking  a  systematic  approach  to  the  problem  of  movement 
artifact  detection  in  functional  echo-planar  imaging. 

Conclusion 

Movement  of  the  imaged  subject  is  a  serious  problem  in 
the  application  of  functional  MR  imaging.  An  attempt  was  made 
in  this  study  to  examine  the  ability  to  detect  movement  by 
analysis  of  echo-planar  images  themselves.  The  algorithms 
developed  and  tested  in  this  study  were  able  to  offer  a 
limited  solution  to  the  problem  of  movement  detection. 
Movement  in  the  plane  of  imaging  could  be  detected  and 
quantified.  Movement  out  of  the  plane  of  imaging,  on  the  other 
hand,  may  at  best  be  detected  only,  with  no  ability  to 
quantitate  it.  The  inability  to  quantitate  movement  out  of  the 
imaging  plane  seriously  limits  the  ability  to  correct  for 
detected  movement  artifact.  As  such,  other  methods  of  dealing 
with  the  problem,  such  as  the  use  of  external  fiduciary 
markers,  more  comprehensive  3-dimensional  image  sets,  or 
better  patient  restraint,  will  have  to  be  considered. 
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Flexure  20 

Proposed  flow-sheet  for  detection  of  movement.  The  center  of 
mass  technique  is  applied  first  to  detect  whether  any 
movement  has  occurred  or  not.  If  movement  is  detected,  the 
area  algorithm  determines  the  nature  of  the  movement  - 
either  in  the  plane  of  imaging  or  outside  of  the  plane  of 
imaging.  Finally,  if  the  movement  is  determined  to  be  in  the 
plane  of  imaging,  the  variance  of  ratios  technique 
determines  the  direction  of  the  movement. 
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APPENDIX 


Image  analysis  program  as  written  by  author  using  the  Zortech© 
C++  compiler  on  an  IBM  PC. 


^include  <  string. h> 

^include  <stream.hpp> 
^include  <conio.h> 

^include  <stdlib.h> 

^include  <math.h> 

^include  <fg.h> 

^include  <sound.h> 

#define  TRUE  1 
^define  FALSE  0 

^define  vx  256 
#define  vy  256 
^define  colormax  256 

^define  sobelx  0 
^define  sobel_y  1 
^define  prewitt_x  2 
^define  prewitt_y  3 
^define  laplace  l  4 
^define  laplace_2  5 
^define  low_pass  6 
^define  high_pass  7 

struct  slice_uc  {  short  int  *im; 

}; 

struct  slice  c  {  fg_color_t  *im; 
" }; 


char  arrow_matrix[15]  = 

{  0x80, 

OxcO, 
OxeO, 
OxfO, 
0xf8, 
Oxfc, 
Oxfe, 
Oxff, 
Oxfc, 
0xf8, 
0xf8, 
Oxcc, 
0x8c, 
0x06, 
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0x06, 

}; 

fgmsmcursort  arrow  —  {arrowmatrix,  {0,0,7, 14},  0,  14}; 

void  button(char  *label,fg_pbox_t  labelbox.int  x,int  y); 
void  clear_box(int  xl,int  yl,  int  x2,  int  y2); 
void  equate(slice_c  *il,  fg_color_t  *i2,int  countz); 
void  equate2(short  int  *il,  short  int  *i2); 
void  equate3(slice_uc  *il,  short  int  *i2,int  countz); 
void  convert( short  int  *il,fg_color_t  *i2,mt  max); 
void  setzero(short  int  *i2); 

short  int  interp(short  int  *image, short  int  *image2, double  x, double  y, double  z); 

void  clear_line(double  a, double  b, double  c); 

void  read_line(char  buffer[],fg_coord_t  x,fg_coord_t  y); 

void  get_center(slice_uc  *i mage, double  *x, double  *y,int  countz); 

void  get_center_w(short  int  *image, double  *x, double  *y); 

short  int  templ(short  int  *image,int  x,int  y,int  temp); 

void  fill(fg_color_t  *image,int  xa,int  ya,long  int  *tot); 

void  get_area_circ(short  int  *image,long  int  *tot_circ,long  int  *tot_area); 

double  get_ratio(short  int  *work, short  int  *work2); 

void  transform(short  int  *work,  short  int  *work2,  int  horizontal,  int  vertical,  double  rotate); 
double  ab(double  number); 

void  n_to_text(int  1, double  numbr,  char  *buffer2); 

void  subtract(short  int  *work,  short  int*work2,  short  int  *work3); 

void  windo(float  *epi, short  int  *work,int  uplimit.int  lowlimit); 

main() 

{ 


int  vz=  11; 


FILE  *fp; 
char  fname[256]; 

slice_uc  *image; 
image  =  new  slice_uc[vz]; 
slice_c  *image_c; 
image_c  =  new  slice_c[vz]; 


//  ARRAYS  FOR  IMAGE  STORAGE 

//  FG  COLOR  T  ARRAYS  FOR  DISPLAY 


short  int  *work; 

work = new  short  int[vx*vy]; 

short  int  *work2; 

work2=new  short  int[vx*vy]; 

short  int  *work3; 

work3  —  new  short  int[vx*vy]; 

short  int  *work4; 

work4=new  short  int[vx*vy]; 


" 
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float  *epi; 

epi  — new  float[vx*vy]; 
float  *epi2; 

epi2  =  new  float[vx*vy]; 

short  int  *xchg; 

xchg  =  new  short  int[vx*vy]; 

short  int  *xchg2; 

xchg2  =  new  short  int[vx*vy]; 

fg_color_t  *temp; 
temp  =  new  fg_color  t[vx*vy]; 
fg_color_t  *temp2; 
temp2=new  fg_color_t[vx*vy]; 


short  int  *shade; 

shade = new  short  int[vx*vy]; 

long  int  *datax; 
datax  =  new  long  int[22]; 
long  int  *datay; 
datay=new  long  int[22]; 

//  START  GRAPHICS 

if  (fg_init_VESA50  =  =  FGNULL) 

{  cout  <  <  "Can’t  start  flash  graphics. 

fgtermO; 

exit(l); 

}; 


//  assign  global  variables 

char  *load_file; 

load  file  =  new  char[30]; 

imsigned  int  status,  las t_ status  =  Oxffff,  quit =0x0000; 
int  coun tx  ==  0 ,  county  -  0 ,  coimtz = ( v z)  /2 ; 
int  color_num,global  counter; 

int  min, max, max2,loc_max, shift, x_tot,y_tot,noise_md=l, hi  =  3000, lo=0; 
int  x  size  =  3 ,  y  size~3,  rsize=3; 

long  int  totjpixels,tot_pixels2,tot_area,tot_circ,tot_area2,tot_circ2; 

double  xO = double(vx) /2 . 0 , yO  =  double( vy ) /2 . 0 , zO  =  double(vz)/2. 0 ; 

double  a  =  80.0,b=60.0,c  =  30.0; 

double  px=0.0,py  =  0.0,pz=0.0,pd=0.0; 

double  theta_z,  threshold =0.5  ,cm  x , cm_y , cm2_x , cm2_y ,  answer; 

fg  coord  t  x,y; 

fg_box_t  ung  scra,img  scm_2 , load_scm , img  box , mo ve , move2 , load , load2 , sve , intrp , qt , save_box , 
fg  box  t  c_of_m, up, dn, It, rt, noise, edge,edge2, area, sobx,soby,prwx,prwy,lapl,lap2,high,low,rst,stl,st2; 
fg_box_t  box  l, box  2,box  3,nxt,prv,acg,ges,vm, ratio, box, trsh, sub, wind; 
fg_box_t  c  high  up , c_high_dwn, c_lo_up ,  c_lo  dwn,  reset ,  show_2 ,  task , v tmp , cpy ; 
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fghandlet  save_handle; 


//  SET  PALLETTE 

for  (color_num=0;  colornum  <  colormax;  color_num++) 

{ 

fg_setPa^ette(c°l°r_nurn » colornum ,  colornum ,  colornum) ; 

r 

fg_setc°lornum(FG_HIGHLIGHT,255); 
fg_setcolomum(F  GWHITE , 255) ; 

fg_setcolomum(F  G_BL  ACK ,  0) ; 

//  SCREEN  DISPLAY 

button("QUIT",qt,215,10); 
button(  "RESET " , rst ,  150, 150); 
button("MOVE",move,  150, 130); 
button("MOVE2",move2, 150, 1 10); 
button("ELLIPSE",st2, 150,90); 
button("C  of  M",c_of_m,260,150); 
button( "  RATIO " ,  ratio , 260 , 1 3 0) ; 
button  ("FILL",area,260,110); 
button  ("A/C  ",acg,260,90); 
button( "  LOAD " ,  load  ,20,150); 
button("NEXT",nxt,20,130); 
button("PREV",prv,20,110); 
button( " — >  " ,  cpy , 20 , 90) ; 
button("INTERP",intrp,20,70); 
button("LOAD2",load2,70,150); 
button("H  UP",c_high_up,70,130); 
button("H  DN",c_high_dwn,70,110); 
button("L  UP",c_lo_up,70,90); 
button("L  DN",c_lo_dwn,70,70); 
button("RESET",reset,70,50); 
button( "  S  A  VE " ,  sve, 20 , 50) ; 
button( "  WIND " ,  wind, 20, 30) ; 
button("SHW2",show_2,70,30); 
button(  "Noise", noise, 215, 90); 
button("THRESH",trsh,215,70); 
button( "  SUBTRACT " ,  sub  ,215,50); 
button("Lapl",lapl,215,150); 
button("Lap2",lap2,215,130); 
button( "  Low ",  low ,  2 1 5 , 1 1 0) ; 

img_box[FG_Xl]  -  5; 
img_box[FG_Y  1]  =  190; 
img_box[FG_X2]  =  635; 
img_box[FG_Y2]  =  475; 

img_scm[FG  Xl]  =  10;  img_scm_2[FG  X  1  ]  =  375; 
img_scm[FG_Y  1]  =  215;  img_scm_2[FG_Y  1  ]  =  2 1 5 ; 
img_scm[FG_X2]  =  265;  img_scm_2[FG_X2]  =  630; 
img_scm[FG_Y2]  =  470;  img_scm_2[FG_Y2]  =  470; 


. 

' 
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box_l[FG_Xl]  =  5;  box_2[FG_Xl]  =  140;  box_3[FG_Xl]  =  340; 
box_l[FG_Yl]  =  5;  box_2[FG_Yl]  =  30;  box_3[FG_Yl]  =  5; 
box_l[FG_X2]  =  130;  box_2[FG_X2]  =  330;  box_3[FG_X2]  =  635; 
box_l[FG_Y2]  =  185;  box_2[FG_Y2]  =  185;  box_3[FG_Y2]  =  185; 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,img_box, fg.displaybox); 
fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,box_l,  fg.displaybox); 
fg_drawbox(FG_WHITE,FG_MODE_SET,~0,FG_LINE_SOLID,box_2, fg.displaybox); 
fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,box_3,fg  displaybox); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,350, 150,  "AREA  1 : "  .fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROT0,350,130,"AREA2:",fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,350, 1 10,  "CofMl  X: " .fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROT0,350,90,"CofM2Y:", fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~ O.FG  ROTO, 350, 70, "CofM  1 X:", fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~ 0, FG  ROTO, 350, 50, "CofM2Y:", fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROT0,350,30,"VARI  AN:  ".fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~ 0,FG_ROT0,5 10, 150,  "X  Shift: " , fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,510, 130, "  Y  Shift: "  .fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  -  0,FG_ROT0,510, 1 10,  "Z  Rot  : " .fg.displaybox); 
clear_box(420,9,509,171); 

//  MAIN  MENU  LOOP 
if  (fg.msm) 

{ 

fg_msm_setcurpos(20 1,110); 
fg_msm_showcursor() ; 
fgflushO; 
do 

{  status  =  fg_msm_getstatus(&x ,  &y ) ; 
if  (status!  =  last  status) 

{ 

if  (status  &  FG  MSM  LEFT) 

{ 

if  (fg_pt_inbox(rst,x,y)) 

{  srand(15); 

global_counter = 0 ; 

for  (int  county  =  0;county  <vy;county+  +) 
for  (int  countx=0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <  8)  +  countx ; 
shade[index]  =  short  int  (rand()/4096  + 1); 

if  ((countx-x0)*(countx-x0)/(a*a)  +  (county-y0)*(county-y0)/(b*b)-1.0 

0.01) 

{  work[index]  =  short  int(FG_WHITE/shade[index]); 
work2[index]  =  short  int(FG_WHITE/shade[index]);} 
else 

{  work[index]  =  FG_BLACK; 
work2[index]  =  FG  BLACK;  } 

} 

max2— FG_WHITE  + 1; 
loc_max  =  max2; 
x_tot=0;y_tot=0; 


■ 
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convert(  work ,  temp ,  locmax) ; 
convert(  work2 ,  temp2 ,  locmax) ; 
fg_wri  tebox(img_scm ,  temp) ; 
fg_writebox(imgscm_2 ,  temp) ; 

fg_puts(FG_WHITE,FG_MODE_SET,  —0, FG_ROTO, 450, 191,  "ORIGINAL",  fg.displaybox); 

} 


if  (fg_pt_inbox(st2,x,y)) 

{ 

char  a_chr[  3  ] ,  b  chr [  3  ] ,  c  chr [  3  ] ,  *dummy ; 
int  base=  10; 

save_box[FG_Xl]  =  140;load_scm[FG_Xl]=  140; 
save_box[FG_Yl]  =  120;load_scm[FG_Yl]=  120; 
save_box[FG_X2]  =  480;load_scm[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;load_scm[FG_Y2]  =  240; 
save_handle  =  fgsave(savebox); 
fg_fillbox(FG_BLACK,FG_MODE_SET,  —  O.loadscm); 
fg_drawbox(FG_WHITE,FG_MODE_SET,  —  0,FG_LINE_SOLID,load_scm,  fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~  O.FG  ROTO,  150,220,  "X*X/A*A  +  Y*Y/B*B  + 
Z*Z/C*C",fg.  display  box) ; 

fg _puts(FG_WHITE,FG_MODE_SET,  ~  O.FG  ROTO,  150, 190,  "ENTER  A  " .fg.displaybox); 
read_lme(a_chr,320, 190); 

fg_puts(FG_WHITE,FG_MODE_SET, -O.FG  ROTO,  150, 170, "ENTER B 
".fg.displaybox); 
read_line(b_chr,320, 170); 

fg_puts(FG_WHITE,FG_MODE_SET,  -  O.FG  ROTO,  150, 150,  "ENTER  C  ".fg.displaybox); 
read_line(c_chr,320, 150); 

a  =  strtol(a_chr ,  &dummy ,  base) ; 
b  =  strtol(b_chr,  &dummy  .base) ; 

c  =  strtol(c_chr,&dummy  .base); 
fgrestore(savehandle) ; 

}" 


if  (fg_pt_inbox(move,x,y)) 

{  int  dec, sign.horiz, vert; 
double  cos  thz = cos(theta_z) ; 
double  sin  thz = sin( theta  z) ; 
double  new_x,new_z,newj,rot; 
char  c_horiz[3],c_vert[3],c_rot[3],*dummy; 
int  base=  10; 

save_box[FG_Xl]  =  140;load_scm[FG_Xl]  =  140; 
save_box[FG_Yl]  =  120;load_scm[FG_Yl]=  120; 
save_box[FG_X2]  =  480;load_scm[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;load_scm[FG_Y2]  =  240; 
savehandle  =  fg_save(save_box); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  -  O.load  scm); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  —0,FG_LINE_SOLID,load_scrn, fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  —  O.FGROTO,  150,220,  "HORIZONTALLY", 
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fg.displaybox); 

read_line(c_horiz,370,220); 

fg_puts(FG_WHITE ,  FG_MODE_SET  ,~0,  FG_ROTO  ,150,200,"  VERTICALLY " , 
fg.displaybox); 
read_line(c_vert ,  3 70 , 200) ; 

fg_puts(FG_WHITE,FG_MODE_SET,  —  0,FG_ROT0, 150, 180, "ROTATE IN  PLANE", 
fg.displaybox); 
read_line(c_rot,370, 180); 
fg_rest°re(  save_handle) ; 


horiz = strtol(c_horiz,  &dummy ,  base) ; 
vert  =  strtol(c_vert ,  &dummy ,  base) ; 
rot  =  0 . 0 1 745329  3  *strtol(c_rot,  &dummy ,  base) ; 
equate2(work2,work3); 
transform(work3,work2,  horiz,  vert,  rot); 
convert(work2 ,  temp2 ,  loc_max ) ; 
fg_writebox(img_scm_2,temp2); 

} 


if  (fg_pt_inbox(move2,x,y)) 

{  srand(16); 

double  new_x,new_y,new_z; 
int  z_perp  =  0,x_rot=0,y_rot  =  0,base=  10; 
char  c_z[3],c_rx[3],c_ry[3],*dummy; 
save_box[FG_Xl]  =  140;load_scrn[FG_Xl]  =  140; 

save_box[FG_Yl]  =  120;load_scm[FG_Yl]  =  120; 
save_box[FG_X2] =480;load_scm[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;load_scra[FG_Y2]  =  240; 
savehandle  =  fg_save(save_box); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~  0,load_scm); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,load_scrn, fg.displaybox); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0, 150,220,  "PERPEDICULAR", 
fg.displaybox); 
read_line(c_z, 3 70 , 220) ; 

fg_puts(F G_ WHITE , F G_M ODE_S ET ,  ~ O.FG  ROTO,  150,200,  "AROUND  X  " , 
fg.displaybox); 
read_line(c_rx  ,370,200); 

fg j5uts(FG_WHITE,FG_MODE_SET,  ~ 0,FG_ROT0, 150,180, "  AROUND  Y" , 
fg.displaybox); 
read_line(c_ry ,  370 , 1 80) ; 
fg_restore(save_handle) ; 
z_perp  =  strtol(c_z,&dummy,base); 
x_rot  =  strtol(c_rx,&dummy,base); 
y_rot  =  strtol(cry ,  &dummy  ,base) ; 
double  costhx  =  cos(0.017453293*x_rot); 
double  sinthx  =  sin(0.017453293*x_rot); 
double  costhy  =  cos(0.017453293*y_rot); 
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double  sinthy  =  sin(0.017453293*y_rot); 
setzero(work2); 

for  (county  =  0;  county  <  vy;  county  +  +) 
for  (countx  =  0;countx<vx;countx  +  +) 

{  int  index  =  (county  <  <  8)  +  countx; 
shade[index]  =  short  int  (rand()/4096  +  1); 
if  (z_perp  !  =  0) 

{  new_x  =  countx-128; 
new_y  =  county- 128; 
new_z=z_perp;} 
if  (x  rot  !  =  0) 

{  new  x  =  countx-128; 

new_y  =  (county- 128)  *cos  thx ; 
new_z  =  (county- 1 28)  *sinthx ; } 
if  (y_rot  !  =  0) 

{  new_x  =  (countx-128)*costhy; 
new_z=(countx-128)*sinthy; 
new_y  =  county- 12  8 ; } 

if  ((new_x*new_x)/(a*a)  +  (new_y*new_y)/(b*b)  +  (new_z*new_z)/(c*c)  -1.0  <  0.01) 
{  work2[index]  =  short  int(FG_WHITE/shade[index]);} 
else 

{  work2[index]  =  FG  BLACK;  }; 

} 

convert(  work2 ,  temp2 ,  locmax) ; 
fg_writebox(img_scm_2 ,  temp2) ; 

} 

if  (fg_pt_inbox(nxt,x,y)) 

{  countz--; 

if  (countz<  1)  {countz=l;  sound_beep(  16000);}; 
if  (countz >(vz-l))  {countz=vz-l;  sound_beep(  16000);}; 
clear_box(  11,216,264,469); 

equate  3  ( i  mage ,  work ,  countz) ;  ’ 
convert(work ,  temp ,  max ) ; 
fg_writebox(imgscm,temp); 
locniax^max; 

} 

if  (fg_pt_inbox(prv,x,y)) 

{  countz  +  +; 

if  (countz  <  1)  {countz  -  1;  sound_beep(  16000);}; 

if  (countz >(vz-l))  {countz =vz-l;  sound_beep(  16000);}; 

clear_box(l  1,216,264,469); 

equates  (i  mage ,  work ,  countz) ; 

convert(work,  temp,  max); 

fg_wntebox(irngscm ,  image_c[countz] .  im); 

locmax  =  max; 

} 


if  (fgjpt_inbox(load2,x,y)) 
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{ 

img_scm[FG_Xl]  =  10;  img_scm_2[FG_Xl]  =  375; 
img_scm[FG_Y  1]  =  215;  img_scm_2[FG_Y  1]  =  2 15; 
img_scm[FG_X2]  =  265 ;  img_scm_2[FG_X2]  =  630; 
img_scm[FG_Y2]  =  470;  img_scm_2[FG_Y2]  =  470;  char  *load_file_2; 
load_file_2  =  new  char[30]; 
load_file_2  =  "  "; 

char  im_no[3],c_max[3],c_min[3],c_hi[7],c_lo[7],*dummy; 
int  base=  10; 

save_box[FG_Xl]  =  140;load_scm[FG_Xl]  =  140; 
save_box[FG_Yl]=  120;load_scm[FG_Yl]=  120; 
save_box[FG_X2]  =  500;load_scm[FG_X2]  =  500; 
save_box[FG_Y2]  =  240;load_scm[FG_Y2]  =  240; 
savehandle  =  fg_save(save_box); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~  O.loadscm); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,load_scrn,fg.displaybox); 


fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0, 150,200, "ENTERFILE 
NAME",fg.displaybox); 
read_line(load_file,300,200); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0, 150, 180,  "ENTER  FILE  NAME 
2",fg.displaybox); 
read_line(load_file_2,300, 1 80); 
sprintf(fname,"  %s",load_file); 
if  ((fp  =  fopen(fname,"rb"))=  =NULL) 

{  fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROT0,300, 140,  "ERROR  CAN’T  OPEN 
FILE " ,  fg .  displaybox) ; 

} 

else 

{ 

fread(epi ,  sizeo  f(  float) ,  vx*vy ,  fp) ; 

} 

fclose(fp); 

if  (load_file_2  !  =  "  ") 

{  sprintf(fname,"%s",load_file_2); 
if  ((fp  =  fopen(fname,"rb"))=  =NULL) 

{  fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,300, 140,  "ERROR  CAN’T  OPEN 

FILE " ,  fg .  displaybox) ; 

} 

else 

{ 

fread(epi2 ,  sizeo  f(  float)  ,vx*vy ,  fp) ; 

} 

fclose(fp); 

}; 


fg_restore(save_handle) ; 


//  CONVERT  TO  REAL  FLOATING  POINT  -  From  Tony  Adrignollo’s  MR  Program  written 


' 
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for  the  VA  Neuropsychology  Laboratory 

short  mt  tmp; 
xchg  =  (short  mt  *)  epi; 
for  (int  i  —  0;i  <  vx*vy-l  ;i++) 

{  tmp  =  *xchg; 

*xchg  =  *(xchg  + 1); 
if  (tmp  j  |  *xchg) 

*(xchg  +  l)  =  tmp-256; 
else  *(xchg+l)  =  tmp; 
xchg  +  =2; 

} 

short  int  tmp2; 
xchg2  =  (short  int  *)  epi2; 
for  (i  =  0;i  <  vx*vy-l  ;i++) 

{  tmp2  =  *xchg2; 

*xchg2  =  *(xchg2  + 1); 
if  (tmp2  |  |  ,(‘xchg2) 

*(xchg2  +  l)  =  tmp2-256; 
else  *(xchg2  +  l)  =  tmp2; 
xchg2  +  =2; 

} 

//  END  CONVERTION  ROUTINE 
setzero(work4); 

for  (int  county  =  0;county  <vy;county  +  +) 
for  (int  countx  —  0; count x  <  vx;countx  +  +) 

{  int  index  =  (county  <  <  8)  +  countx; 

int  index2=  ((vy-coimty-l)<  <8)  +  countx; 
work[index2]  =  epi[index]; 
work4[index2]  =  epi2[index]; 

} 

clear_box(30, 191, 190,212); 
n_t°  _text(5  ,hi ,  chi) ; 
n_t°_text(5 ,  lo ,  c_lo) ; 

fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,55,192,c_hi,fg.displaybox); 
fgjputs(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,130,192,c_lo,fg.displaybox); 

locrnax  =  9999 ; 
windo(epi  .work, hi, lo); 

convert(work,  temp ,  locmax) ; 
fg_writebox(img_scra,  temp) ; 

windo(epi2 ,  work4 ,  hi ,  lo) ; 
convert(work4 ,  temp2 ,  loc_max) ; 
fg_writebox(img_scm_2 ,  temp2) ; 


} 


if  (fgjpt_inbox(reset,x,y)) 
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{  windo(epi,work,hi,lo); 
loc  max  =  9999; 
convert(work ,  temp ,  locmax) ; 
fg_writebox(img_scra ,  temp) ; 

equate2(work,  work2) ; 
convert(work2 ,  temp  2 ,  locmax) ; 
fg_writebox(img_scm_2,temp2); 

} 

if  (fg_pt_mbox(c_high_up,x,y)  j  j  fg_pt_inbox(c_high_dwn,x,y)  {  j  fg_pt_inbox(c_lo_up,x,y) 
fg_pt_inbox(c_lo_dwn ,  x ,  y )) 

{char  c_hi[7],c_lo[7]; 
if  (fg_pt_inbox(c_high_up,x,y)) 
hi  =  hi  + 100; 

if  (fg_pt_inbox(c_high_dwn,x,y)) 
hi  =  hi-100; 

if  (fg_pt_mbox(c_lo_up,x,y)) 
lo  =  lo  + 100; 

if  (fg_pt_inbox(c_lo_dwn,x,y)) 
lo  =  lo-100; 

clear_box(30, 191, 190,212); 
n_to_text(5,hi,c_hi); 
n_to_text(5 ,  lo ,  c_lo) ; 

fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROTO,55,192,c_hi,fg.displaybox); 
fg_puts(F  G_WHITE ,  F  G_M  ODE_SET,~0,F  GROTO ,  1 30 , 1 92 ,  c_lo ,  fg .  display  box) ; 

windo(epi ,  work  ,hi ,  lo) ; 
loc_max  =  9999; 

convert(  work ,  temp ,  locmax) ; 
fg_writebox(img_scm,  temp); 


} 


if  (fg_pt_inbox(show_2,x,y)) 

{  windo(epi2,work4,hi,lo); 
loc_max  =  9999; 
convert(work4 ,  temp2 ,  locmax) ; 
fg_writebox(img_scm_2 ,  temp2); 
equate2(work4 ,  work2) ; 

} 


//■ 


if  (fg  pt  mbox(sve,x,y)) 

{ 
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char  *save_file; 
save  file  =  new  char[30]; 
char  im_no[3],frmt[l],*dummy; 
int  base  =  10; 

save_box[FG_Xl]  =  140;load_scm[FG_Xl]=  140; 
save_box[FG_Yl]  =  120;load_scm[FG_Yl]  =  120; 
save_box[FG_X2]  =  480;load_scm[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;  load_scm[FG_Y2]  =  240; 
save_handle  =  fgsave(savebox); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~ O.load  scra); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,load_scm,fg.displaybox); 


fg_puts(FG_WHITE,FG_MODE_SET,  ~  O.FGROTO,  150,200,  "ENTER  FILE 
NAME",fg.displaybox); 

read_line(save_file,300,200); 

sprintf(fhame , "  %  s .  00 1 " ,  save_file) ; 
if  ((fp  =  fopen(  fname, "  wb " ))  =  =  NULL) 

fg_puts(FG_WHITE,FG_MODE_SET,  ~  0, FG  ROTO, 300, 140,  "ERROR  Cannot  Open 

File ",  fg .  display  box) ; 

else 

fwrite(work2,sizeof( short  int),vx*vy,fp); 
fclose(fp); 

fgrestorefsayehandle) ; 

} 

if  (fg_pt_inbox(load,x,y)) 

{ 


char  im_no[3],frmt[l],*dummy; 
int  base=  10; 

save_box[FG_Xl]=  140;load_scm[FG_Xl]=  140; 
save_box[FG_Yl]  =  120;load_scm[FG_Yl]  =  120; 
save_box[FG_X2]  =  480;load_scrn[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;load_scrn[FG_Y2]  =  240; 
savehandle  =  fgsave(savebox); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~  0,load_scrn); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,load_scra,fg.displaybox); 
fg_puts(F G^WHITE , F G_M ODE  SET ,  ~  0,FG_ROT0, 150,200,  "ENTER  FILE 
NAME",fg.displaybox); 

read_line(load_file,  3 00 , 200) ; 

fg_puts(F G_WHITE , F G_M ODE_SET ,  ~  O.FG  ROTO,  150,1 80,  "NUMBER  OF 
IMAGES  ",fg.displaybox); 

read_line(ixn_no,300, 180); 

vz = int(strtol(im_no,  &dummy  ,base)) ; 

for  (countz=l;countz<vz;countz+  +) 

{  sprintf(fname, " %s. %03d",load_file,countz); 
if  ((fp  =  fopen(fhame,"rb"))=  =NULL) 


* 
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{  fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,300, 140,  "ERROR  Can’t  Open 

File",fg.displaybox); 

exit(0);} 

else 

{  image[countz].im=new  short  int[vx*vy]; 
fread(work,  si  zeof(short  int) , vx*vy ,  fp) ; 
for  (county  =  0;county <vy;county  +  +) 
for  (countx  =  0;countx<  vx;countx+  +) 

{  int  index  =  (county  <  <8)  +  countx; 
int  index2  =  ((vy-county-l)<  <8)  +  countx; 
image[countz] .  im[index2]  =  work[index] ; 

}; 

fclose(fp);} 

} 

//  GET  MAX/MIN  FOR  SCALING 


max  =  0,  min  — 9999; 


for  (countz=  l;countz<  vz;countz+ +) 
for  (county  =  0;county  <vy;county  +  +) 
for  (countx  =  0;countx<vx;countx  +  +) 

{int  index  =  (county  <  <  8)  +  countx; 

if  (image[countz].im[index]<min)  min  =  image[countz].im[index]; 
if  (image[countz].im[index]  >  max)  max  =  image[countz].im[index]; 


}; 

locmax  =  max ; 

//  RE-SCALE  &  FLIP  IMAGE  TO  MAKE  ORIGIN  AT  BOTTOM  LEFT  -  MAKE  "IMAGE  C" 

for  (countz=  l;countz<  vz;countz+  +) 

{image_c[countz] ,im= new  fg_color_t[vx *vy ] ; 
for  (county =0;county  <vy;county+  +) 
for  (countx  =  0;countx  <  vx;countx+  +) 

{  int  index  =  (county  <  <  8)  +  countx; 

image_c[countz].im[ index]  =  (fg  color  t)  (image[countz] . im( index]*color_max/max); 

} 

} 

fgrestore(savehandle) ; 

countz=(vz)/2; 

equate3  (i  mage,  work ,  countz) ; 

convert(  work ,  temp ,  max) ; 

fg_writebox(img_scm,image_c[countz] .  im); 
fgwritebox(img_scm_2,temp); 

} 


if  (fgj>t_inbox(noise,x,y)) 

{  fg  color  t  i  =  1 ; 
int  j  =  2000; 

for  (int  county  =  0; county  <vy;county+  +) 
for  (int  countx =0;countx<  vx; countx  +  +) 
{  srand(rand()-county  +  noise_md); 
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int  index  =  (county  <  <8)  +  countx; 
if  (rand()<j)  work2[index]  =  loc_max-l; 

}; 

convert(work2 ,  temp2 ,  locmax) ; 
fg_writebox(img_scm_2 ,  temp2) ; 
noise_md+  + ; 

} 

if  (fg_pt_inbox(intrp,x,y)) 

{  char  c  depthf  3  ] ,  *dummy ; 
double  depth  =  0, base  =  10; 
save_box[FG_Xl]  =  140;load_scm[FG_Xl]  =  140; 
save_box[FG_Yl]  =  120;load_scm[FG_Yl]  =  120; 
save_box[FG_X2]  =  500;load_scm[FG_X2]  =  500; 
save_box[FG_Y2]  =  240;load_scm[FG_Y2]  =  240; 
savehandle  =  fg_save(save_box); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~  O.loadscm); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,load_scm,fg.displaybox); 


fg_puts(FG_WHITE , FG_M ODE  SET ,  ~  0,FG_ROT0, 150,200,  "DIST  FROM  IMG 
l",fg.  display  box); 

read_line(c_depth,300,200); 

depth=double(strtol(c_depth,&dummy,base)),,,0. 1428571; 
fgrestore(savehandle) ; 


for  (int  county  =  0;county  <vy; county  +  +) 
for  (int  countx=0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <8)  +  countx; 

work2[index]  =  interp(work,  work4 ,  double(countx) ,  double(county) ,  depth) ; } 
convert(work2 ,  temp2 ,  locmax) ; 
fg_writebox(img_scm_2 ,  temp2) ; 

} 


if  (fg_pt_inbox(c_of_m,x,y)) 

{  char  c_cmxl[5],c_cmyl[5],c_cmx2[5],c_cmy2[5],c_xdiffI5],c_ydifft5]; 
get_center_w(  work ,  &cm_x ,  &cm_y ) ; 
get_center_w(work2  ,&cm2_x ,  &cm2_y) ; 
int  dec, sign; 

clear_box(420 , 49 , 509 , 1 3 1 ) ; 
clear_box(574, 109,634, 171); 

n_to_text(3  ,cmx,c_cmx  1 )  ;n_to_text(3 ,  cm_y ,  c_cmy  1 ) ; 
n_to_text(3 ,  cm2_x ,  c_cmx2)  ;n_to_tex  t(  3 ,  cm2_y ,  c_cmy  2) ; 
n_t0_text(3,cm2_x-cm_x,c_xdiff);n_to_text(3,cm2_y-cm_y,c_ydiff); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,425, 1 10,c_cmxl  ,fg.  display  box); 
fg_puts(F  G_WHITE ,  F  G_M  ODE_SE*T ,  ~0,F  GROTO , 425 , 90 ,  c_cmy  1 ,  fg .  di  splay  box) ; 
fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,425,70,c_cmx2,fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ^~0,FG_ROT0,425,50,c_cmy2,fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  —  0,FG_ROT0,575,150,c_xdiff,fg.displaybox); 


' 

* 
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fg_puts(FG  WHITE, FG_MODE_SET,  ~  0,FG_ROT0, 575, 130,c_ydiff,fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,575,110,"  ?",fg.displaybox); 


if  (fg_pt_inbox(trsh,x,y)) 

{ 

char  c_thresh[3],*dummy; 
int  base  =  10; 

save_box[FG_Xl]  =  140;load_scm[FG_Xl]  =  140; 
save_box[FG_Y  1]  =  120;load_scrn[FG_Y  1]  =  120; 
save_box[FG_X2] =480;load_scm[FG_X2]  =  480; 
save_box[FG_Y2]  =  240;load_scm[FG_Y2]  =  240; 
savehandle  =  fg_save(save_box); 
fg_fillbox(FG  BLACK, FG  MODE  SET,  ~ 0,load_scrn); 

fg_drawbox(F  GWHITE ,  F  GMODESET ,  ~  0 ,  FG_LINE_S  OLID ,  load_scm,  fg .  display  box) ; 
fg j3uts(FG_WHITE,FG_MODE_SET,  ~ O.FG  ROTO,  150,220, "NEW 
THRESHOLD",  fg.displaybox); 

read_line(c_thresh,  370 ,220) ; 
fg_restore(save_handle); 

threshold  =  double(strtol(c_thresh ,  &dummy ,  base))/ 100.0; 


} 

if  (  fg_pt_inbox(lapl,x,y)  |  |  fg_pt_inbox(lap2,x,y)  |  |  fg_pt_inbox(low,x,y)) 

{ 

short  int  local_max  =  0,local_max2  =  Q; 
for  (county  =  0;county  <  (vy -2); county  +  +) 
for  (countx  =  0;countx  <  (vx-2);countx+ +) 

{  int  index  =  (county  <  <8)  +  countx; 

if  (fg_pt_inbox(lapl,x,y))  {  work2[index]  =  templ(work2,countx, county, laplace  l); 

workfindex]  =  templ(work,  countx ,  county ,  laplace_  1 ) ; } 
if  (fg_pt_inbox(lap2,x,y))  {  work2[index]  =  templ(work2, countx, county, laplace_2); 

work[index]  =  templ(work,  countx ,  county ,  laplace_2) ; } 
if  (fg_pt_inbox(low,x,y))  {  work2[index]  =  templ(work2, countx, county, low_pass); 

workfindex]  =  templ(work, countx, county, low_pass);} 
if  (local  max  <  work[index])  local_max=work[index]; 
if  (local_max2  <  work2[index])  local_max2  =  work2[index]; 

} 

//  THRESHOLD  VALUES 

for  (county  ==0; county  <  vy; county  +  +) 
for  (countx  =  0; countx  <  vx ; countx -+-  +) 

{  int  index  —  (county  <  <  8)  +  countx; 
if  (work[index]  <  (short  mt)(threshold*local  max) )  work[index]  =  0; 
else  work[index]  =  (short  int)local  max-l; 

if  (work2[index]  <  (short  int)(thre«hold*local_max2) )  work2[index]  =  0; 
else  work2.[iridex]  ~  (short  int)local_max2-l; 

} 

convert(  work,  temp ,  local_max  + 1 ) ; 
convert(work2,temp2,local_max2  + 1); 
fgwritebox(imgscm,  temp) ; 


' 

■ 


Page  93 


fg_writebox(img_scm_2 ,  temp2) ; 


} 

if  (fgjpt_inbox(sub,x,y)) 

{  long  int  tot_ar=0; 
char  c_ar[7]; 

subtract(work2,work,work3); 

equate2(work3  ,work2); 
convert(work2 ,  temp2 ,  locinax); 
fg_writebox(img_scm_2,temp2); 

for  (county  -  0 ;  county  <  vy ;  county  +  +) 
for  (countx  =  0;countx<vx;countx  +  +) 

{  int  index  =  (county  <  <8)  +  countx; 
if  (work2[index]  !=  0)  tot_ar+  +; 

} 

n_to_text(5 ,  tot_ar  ,c_ar) ; 
clear_box(420, 129,509, 171); 

fg_Puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,430,130,c_ar,fg.displaybox); 

} 

if  (fg_pt_inbox(area,x,y)) 

{  char  c_totl[7],c_tot2[7]; 
totpixels  —  0;  tot_pixels2  =  0; 
fill(temp,  128, 128,&tot_pixels); 
fill(temp2, 128, 128,&totjpixels2); 
fg_writebox(img_scm_2 ,  temp2) ; 

fg_writebox(img_scrn ,  temp) ; 
int  dec, sign; 

n_to_text(5 ,  totjhxels ,  c_tot  1 ) ;  n_to_text(5 ,  tot_pixels2 ,  c_tot2) ; 
clear_box(420, 129,509, 171); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~  0,FG_ROT0,430, 150,c_totl  ,fg.displaybox); 
fg_puts(FG_WHITE,FG_MODE_SET,  ~0,FG_ROT0,430,130,c_tot2,fg.displaybox); 

}; 


if  (fg_pt_inbox(ratio,x,y)) 

{  answer =0.0; 
double  min ;  min  =  99999 . 0; 

char  c_hbas[3],c  vbas[3],c_rbas[3],c  hp[3],c  vp[3],c_rp[3],+dunimy; 
char  c_hor[5],c  ver[5],c  rot[5],c_min[7]; 
int  base  =10; 


save_box[FG_Xl]  =  140;load_scm[FG_Xl]=  140; 
save_box[FG_Y  1]  =  120;loadscm[FGY  1]  =  120; 
save_box[FG_X2]=480;load_scm[FG_X2]  =480; 
save_box[FG  Y2]  =  240;load_scrn[FG_Y2]  =  240; 
save_handle  =  fg_save(save_box); 

fg_fillbox(FG_BLACK,FG_MODE_SET,  ~0,load_scm); 

fg_drawbox(FG_WHITE,FG_MODE_SET ,  —  0 ,  F  GLINESOLID ,  loadscm,  fg .  display  box) ; 


.. 
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fg_puts(FG_WHITE,FGJMODE_SET,  ~  O.FGROTO,  150,220,  "BEST  X 
GUESS", fg.displaybox); 

read_line(c_hbas  ,370,220); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~  O.FGROTO,  150,200,  "BEST  Y 
GUESS", fg.displaybox); 

read_line(c_vbas ,  370, 200) ; 

fg_puts(FG_  WHITE, FGMODESET,  ~  O.FGROTO,  150, 1 80,  "BEST  Z 
GUESS",  fg.displaybox); 

read_l  ine(c_rbas  ,370,180); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~ O.FG  ROTO,  150, 160, "X points" .fg.displaybox); 
read_l  ine(c_hp  ,370,160); 

fg_puts(FG_WHITE,FG_MODE_SET,~0,FG_ROTO,  150, 140,  "Y  points",  fg.displaybox); 
read_line(c_vp,370, 140); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~  O.FG  ROTO,  150, 120,  "Rotations" .fg.displaybox); 
read_line(c_rp,370, 120); 
fg_restore(save_handle) ; 

int  x_size  =  strtol(c_hp,&dummy,base); 
int  y  size  =  strtol(c_vp,&dummy,base); 
int  r_size  =  strtol(c_rp,&dummy,base); 

int  base_h  =  strtol(c_hbas,&dummy,base); 
int  base_v  =  strtol(c_vbas ,  &dummy ,  base) ; 
int  base_r  =  strtol(c_rbas,&dummy,base); 

double  *square; 

square = new  double[x_size*y_size*r_size] ; 

int  xs2  =  (x_size-l)/2,  ys2  =  (y_size-l)/2,rs2  =  (r_size-l)/2; 

equate2(work2 ,  work3 ) ; 

int  hor,ver,rotat,dec,sign;hor=0;ver=0;rotat  =  0; 

for  (int  h=0;h  <  x_size;h+  +) 
for  (int  v  =  0;v  <  y_size;v+  +) 
for  (int  r=0;r  <  r_size;r+  +) 

{  int  index  =  (h*r_size*y  size) +  (v*r_size)  +  r; 
transform(work ,  work2 ,  baseh  +  h-x  s2 ,  base_v  +  v-y  s2 ,  baser  +  (r-rs2)  *0.017453293); 
square[index]  =  get_ratio(work2,work3); 
clear_box(419,9,560,91); 

char  c_count[5]  ;n_to_text(3 ,  double(index) ,  c  count); 

fg_puts(FG_WHITE,FG_MODE_SET,  ~0, FG  ROTO, 420, 10, "COUNTER: ".fg.displaybox); 
fg_puts(F G  WHITE , F G_M ODE  SET ,  ~ 0,FG_ROT0,490, 10,c_count,fg.displaybox); 

}; 


for  (h=0;h<x_size;h+  +) 
for  (v  =  0; v  <  y_size; v  +  + ) 
for  (r=0;r<r_size;r+  +) 

{  int  index  =  (h*r_size*y  size) +  (v*r_size)  +  r; 
if  (square[index]  <  min)  {  min  —  square[index]; 

hor  =  base  h  +  h-xs2 ; ver = base_v  +  v-ys2 ;  rotat  =  base  r  +  r-rs2;  } 
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} 

answer = min; 

clear_box(574, 109,634, 171); 
clear_box(4 19,9,634,51); 

n_to_text(3>hor,c_hor);n_to_text(3,ver,c_ver);n_to_text(3,rotat,c_rot); 

n_to_text(5 ,  min,  cmin) ; 


}; 


if  (fg_pt_inbox(cpy,x,y)) 

{  equate2(work,work2); 

convert(work2 ,  temp2 ,  locmax); 
fg_writebox(img_scm_2 ,  temp2) ; 

}; 


if  (fg_pt_inbox(qt,x,y)) 

{ 

quit  =  1; 
break; 

} 


} 

} 

fg_flush(); 
laststatus = status; 

}  while  (!quit); 


} 

//  END  GRAPHICS 
fgflush; 

fgtermO 
exit(O); 
return  0; 


}; 


void  button(char  *label,fg  box  t  label  box,int  x,int  y) 

{ 

fg_coord_t  xc  =  fg_coord_t(x); 
fg  coord  t  yc  --  fg_coord_t(y); 
size  t  len  =  strlen(label); 

1  abelbox  [FGJX 1  ]  =  xc- 1 ; 
label_box[FG_Y  1]  =  yc; 

label_box  [F  G_X2]  =  xc  +  len*fg_box_width(fg.charbox); 
label_box[FG_Y  2]  =  yc  +  fgboxheightlfg .  charbox); 

fg_drawbox(FG_WHITE,FG_MODE_SET,  ~0,FG_LINE_SOLID,label_box,fg.displaybox); 
fg__puts(F  GWHITE  ,FG_M  ODE_SET ,  ~0,F  GROTO ,  xc ,  yc ,  label ,  fg .  display  box) ; 


■ 
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}; 


void  clear_box(int  xl.int  yl,  int  x2,  int  y2) 

{  fg_coord_t  xlc  =  fg_coord_t(xl); 
fgcoordt  x2c  =  fg_coord_t(x2); 
fgcoordt  ylc  =  fg_coord_t(yl); 
fg  coord  t  y2c  =  fg_coord_t(y2); 
fg  box  t  black  box; 
black_box  [F  G_X  1  ]  =  x  1  c ; 
black_box[FG_X2]  =  x2c; 
blackbox  [FG_Y  1  ]  =  y  1  c ; 
black_box[FG_Y2]  =  y2c; 

fg_fxllbox(FG_BLACK,FG_MODE_SET,  -  0,black_box); 

} 

void  equate(slice_c  *il,  fg_color_t  *i2,int  countz) 

{  for  (int  county  =  0;county  <  vy ;county  +  + ) 

for  (int  countx  =  0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <  8 )  +  countx ; 
i2[index]  =  i  1  [countz] .  im[  index] ; 


void  equate2(short  int  *il,  short  int  *i2) 

{for  (int  county  =  0;county  <vy;county+  +) 

for  (int  countx=0;countx<vx;countx+  +) 
{  int  index  =  (county  <  <  8)  +  countx ; 
i2[index]  =  i  1  [index] ; 

}; 

}; 


void  equate3(slice_uc  *il, short  int  *i2,int  countz) 
{for  (int  county  =  0; county  <vy; county  +  +) 

for  (int  countx  =  0;countx<vx;countx+  +) 
{  int  index  =  (county  <  <  8)  +  countx; 
i2[index]  =  i  1  [countz] .  im[  index] ; 

}; 

}; 


void  setzero(short  int  *i2) 

{  for  (int  county =0;county  <vy;county+  +) 

for  (int  countx =0;countx<  vx;countx+  +) 
{  int  index  =  (county  <  <  8)  +  countx ; 
i2[index]  =  0; 

}; 

}; 


void  read_line(char  buffer[],fg_coord_t  x,fg_coord_t  y) 
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char  character; 

int  i  =  0; 
do 

{  character=getch  (); 

if  (character  !  =  ’\n’  &&  character  !=  ’\r’  &&  character  !  =  ’\b’) 

fg_putc(FG_WHITE,FG_MODE_SET,  -'0,FG_ROT0,x  +  i*fg_box_width(fg.charbox),y,  character,  fg.  display  box); 
if  (character  =  =  ’\r’) 
character  =  ’\n’; 
buffer[i]  =  character; 

+  +  i; 


} 

while  (character  !  =  ’\n’); 
buffer[i-l]  =  ’\0’; 


void  get_center(slice_uc  '•‘image, double  *x, double  *y,int  coimtz) 

{  long  int  x_tot  =  0,y_tot=0,z_tot=0,m_tot  =  0; 

*x  =  0;*y  =  0; 

for  (int  county  =  0;county  <vy;county+  +) 
for  (int  countx  =  0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <  8)  +  countx; 

xtot  +  =  countx  *image[countz] .  im[  index] ; 
y_tot+  =county*image[countz].im[index]; 
m_tot+  =image[countz].im[index]; 

} 

if  (m_tot!  =0) 

{  *x=(x_tot/m_tot); 

*y  =  (y_tot/m_tot) ; 

}  ' 
else 

{  *x  =  0.0; 

*y  =  0.0;} 

} 

void  get_center_w(short  int  *image, double  *x, double  *y) 

{  long  int  x_tot  =  0,y_tot=0,z_tot=0,m_tot=0; 

*x=0;*y=0; 

for  (int  county =0;county  <vy;county+  +) 
for  (int  countx=0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <  8)  +  countx ; 
x_tot+  —  countx  *image[index] ; 
y_tot+  =  county ’l‘image[index]; 
mtot  +  =image[index]; 

} 

if  (m_tot!  =0) 

{  *x  =  double(x_tot/m_tot); 

*y  =  double(y_tot/m_tot); 

} 

else 
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{  *x=0.0; 

*y  =  0.0; 

}; 

} 

void  convert(short  int*  il,fg_color_t  *i2,int  max) 

{ 

for  (int  county  =  0;county  <  vy;county  +  +) 
for  (int  countx  =  0;countx<vx;countx  +  +) 

{  int  index  =  (county  <  <8)  +  countx; 
i2[index]  =  (fg_color_t)  ((il[index]*color_max)/max); 

} 

} 

short  int  interp(short  int  *image, short  int  *image2, double  x, double  y, double  z) 
{  //  As  adapted  from  Coritech’s  interpolation  routine 
int  ret; 
int  ix,iy,iz; 

ix  =  int(x);  iy  =  int(y);  iz=int(z); 

if  (ix<0  ||  ix>255  |  |  iy<0  |  |  iy>255  |  |  iz<0  |  |  iz>2) 
return  0; 

int  dxl,dyl,dzl,dx2,dy2,dz2,index  =  (iy <  <8)  +  ix; 

dxl  =  16-(dx2  =  int((x-ix)*16)); 

dyl  =  16-(dy2  =  int((y-iy)*16)); 

dzl  =  16-(dz2  =  int((z-iz)*16)); 

ret = dx  1  *dy  1  *dz  1  *image[index] ; 

ret  +  =  dx  1  *dy  1  *dz2  *i  mage2  [index  ] ; 

ret+  =  dxl*dy2*dzl*image[index  +  256]; 

ret+  =dxl+dy2*dz2*image2[index  +  256]; 

ret  +  =  dx2*dy  1  *dz  1  *image[index  + 1  ] ; 

ret+  =dx2*dyl*dz2*image2[index  + 1]; 

ret+  =dx2*dy2*dzl*image[index  +  257]; 

ret  +  =  dx2*dy2*dz2*image2[index  +  257] ; 

ret>  >  =  12; 

return  (ret<0)?0:(ret); 


} 

short  int  templ(short  int  %image,int  x,int  y,int  temp) 

{ 

int  xl  =  0,x2  =  0,x3  =  0,yl=0,y2=0,y3  =  0,zl=0,z2  =  0,z3  =  0; 


if  (temp=  =sobel_x) 

{  xl  =  -l  ;  x2=  0;  x3=  1; 
yl  =  -2  ;  y2=  0;  y3=  2; 
zl  =  -1  ;  z2=  0;  z3=  1; 

} 

if  (temp  =  =sobel_y) 

{  xl=  1  ;  x2=  2;  x3=  1; 
yl=  0  ;  y2=  0;  y3=  0; 
zl  =-l  ;  z2  =  -2;  z3  =  -l; 

} 
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if  (temp  =  =prewitt_x) 

{  xl  "  -1  ;  x2  =  0;  x3  =  1; 
yl  =  -1  ;  y2=  0;  y3  =  1; 
zl  =  -1  ;  z2=  0;  z3  =  1; 

} 

if  (temp  =  =prewitt_y) 

{  xl—  1  ;  x2—  1;  x3==  1; 
y  1  —  0  ;  y2=  0;  y3=  0; 
zl  =  -l  ;  z2  =  -l;  z3--l; 

} 

if  (temp=  -laplace_l) 

{  xl=  0  ;  x2  =  -l  ;  x3=  0; 
yl=-l  ;  y2=  4  ;  y3=-l; 
zl=  0  ;  z2  =  -l  ;  z3=  0; 

} 

if  (temp=  =laplace_2) 

{  xl  =  -l  ;  x2  =  -l;  x3  =  -l; 
yl  =  -l  ;  y2=  8;  y3  =  -l; 
zl  =  -l  ;  z2  =  -l;  z3  =  -l; 

} 

if  (temp=  =high_pass) 

{  xl=  0  ;  x2  =  -l;  x3  =  0; 
yl=-l  ;  y2=  4;  y3  =  -l; 
zl -=  0  ;  z2  =  -l;  z3=  0; 

} 

if  (temp  —  =low_pass) 

{  xl  =  1;  x2=  1;  x3=  1; 
yl=  1;  y2=  1;  y3  =  1; 
zl  =  1;  z2=  1;  z3=  1; 

}; 

int  ret=0; 

int  index =  (y  <  <8)+x; 
ret = x  1  *image[index] ; 
ret  +  =  x2*image[index+  1]; 
ret  +  =  x3  *image[index  +  2] ; 
ret  +  =  yl*image[index  +  256]; 
ret  +  =  y2*image[index+257]; 
ret+  =  y3*image[index  +  258]; 
ret  +  =  zl  *image[index  +  5 12]; 
ret  -I-  =  z2*image[index  +  513]; 
ret+  =  z3*image[index  +  514]; 
return  ret; 

} 

void  fill(fg_color_t  *image,int  xa,int  y a, long  int  *tot) 
{  int  xleft.xright; 
int  oldxa  -xa; 
int  index  =  (ya<  <8)  +  xa; 

//  fill  to  left 

while  (imagefindex  j  =  =  0  &&  xa  >  0) 
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{ 

image[index]  =  FGWHITE; 

*tot+  =  l; 
xa--; 

index  =  (ya  <  <  8)  +  xa; 

} 

xleft  =  xa+ 1; 

//  fill  to  right 

xa  =  oldxa+ 1; 

index  =  (ya<  <8)  +  xa; 

while  (image[index]  =  =  0  &&  xa  <  vx) 

{ 

imagefindex]  =  FG_  WHITE; 

*tot+  =  1; 
xa  +  +; 

index  =  (ya  <  <  8)  +  xa; 

} 

xright  =  xa; 

//  Check  above  row  &  recursive  fill  of  new  area 
for  (xa=xleft;xa  <  xright;xa++) 

{  int  index  =  ((ya  + 1)  <  <  8)  +  xa; 
if  (ya<vy  &&  image[index]  =  =0)  fill(image,xa,ya+  l,tot); 

} 

//  Check  below  row  &  recursively  fill  new  area 
for  (xa=xleft;xa  <  xright;xa++) 

{  int  index  =  ((ya-l)<  <8)  +  xa; 

if  (ya>0  &&  image[index]=  =0)  fill(image,xa,ya-l,tot); 

} 

}; 


double  get_ratio(short  int  *work, short  int  ♦world) 

{  int  n, dec, sign, criti; 

double  tot  r,  tot  st  dv ,  mean  r ,  std_dev_r ,crit; 

tot_r=0.0;mean_r=0.0;tot_st_dv  =  0.0;std_dev_r=0.0;crit=0.0;n=0 

double  *ratio; 

ratio = new  doublet vx*vy]; 

for  (int  county  =  0;county  <vy; county  4-  +) 
for  (int  countx  =  0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <8)  +  countx; 
ratio[index]  =  double((work[index]  +  l)/(work2[index]  + 1)); 
n  +  +; 

tot_r+  =ratio[index]; 

} 

mean_r = totr/ double(n) ; 
for  (county  =  0;county  <vy;county+  +) 
for  (countx=0;countx<vx;countx+  +) 

{  int  index  =  (county  <  <8)+countx; 
totstdv  +  =  (ratio[index]-mean_r)*(ratio[index]-mean_r); 

} 
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std_dev_r = sqrt(tot_st_dv/ (double(n)- 1.0)); 
if  (std_dev_r!=0.0) 
crit  =  mean_r/std_dev_r; 
else 

crit  =  0.0; 
delete  ratio; 
return  crit; 

} 

void  transform(short  int  *work,  short  int  *work2,  int  horizontal,  int  vertical,  double  rotate) 
{  double  cos  thz  =  cos( rotate); 
double  sin  thz  —  sin(rotate); 
double  new_x,new_y,new_z; 

for  (int  county  =  0;county  <vy;county  +  +  ) 
for  (int  countx  =  0;countx<vx;countx  +  +  ) 

{  int  index  =  (county  <  <  8)  +  countx; 
work2[index]  =  0; 

}; 


int  roto=  128; 

for  (county =0;county  <vy; county  +  +) 
for  (countx  =  0;countx<vx;countx  +  +) 

{  int  index  =  (county  <  <8)  +  countx; 
newx  =  (countx-roto)*cos_thz-(county-roto)*sin_thz  +  horizontal; 
new_y  =  (countx-roto)*sin_thz+ (county-roto)*cos_thz  +  vertical; 
new_z=0; 

int  index  new  =  int(new_y  +  roto)  *vy  +  int(new_x  +  roto) ; 
if  ((index_new  >  vx*vy)  |  |  (index_new)  <  0) 

{} 

else 

work2[index_new]  =  workfindex] ; 

>; 


double  ab(double  number) 

{ 

if  (number  <  0)  number = -number; 
return  number; 

}; 


void  n_to_text(int  1, double  numbr,  char  *buffer2) 

{ 

char  *buffer; 

//  char  buffer2[l  +  2]; 
int  dec, sign; 

buffer = ecvt(numbr ,  1 ,  &dec ,  &sign) ; 

if  (sign  !=  0)  buffer2[0]  =  ’-’;  else  buffer2[0]  =  ’  ’; 
for  (int  i=0;i<l  +  l;i+  +) 

{if  (i  <  dec)  buffer2[i  +  l]= buffer [i]; 
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if  (i  ==  dec)  buffer2[i  + 1]= 

if  (i  >  dec)  buffer2[i+  l]  =  buffer[i-l]; 


}; 


void  subtract( short  int  *work,  short  int*work2,  short  int  *work3) 

{ 

for  (int  county  =  0;county  <vy;county  +  +  ) 
for  (int  countx  =  0;countx<  vx;countx  +  +) 

{  int  index  =  (county  <  <8)  +  countx; 
work3  [index]  =  work[index]-work2[index]; 

}; 

}; 


void  windo(float  *epi, short  int  *work,int  uplimit.int  lowlimit) 

{ 

for  (int  county  =  0;county  <vy; county  +  +) 

for  (int  countx  =  0;countx<  vx;countx+  +) 

{int  index  =  (county  <  <8)  +  countx; 
int  index2  =  ((vy-county-l)<  <8)  4-  countx; 
if  (epi[index]  <  FG  BLACK) 
work[index2]  =  lowlimit; 
else  if  (epi[index]  >  uplimit) 
work[index2]  =  FGBLACK; 

else  work[index2]  =  ((epi[index]-lowlimit)/(uplimit-lowlimit) *9999) ; 


} 


} 
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